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ABSTRACT 


A  modeling  technique  has  been  developed  to 
account  for  interactions  between  load  current  and  magnetic 
flux  in  an  iron  conductor.  Such  a  conductor  would  be  used 
in  the  active  region  of  a  normally  conducting  homopolar 
machine.  This  approach  has  been  experimentally  verified 
and  its  application  to  a  real  machine  demonstrated. 
Additionally,  measurements  of  the  resistivity  of  steel  under 
the  combined  effects  of  magnetic  field  and  current  have 
been  conducted. 
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BACKGROUND 


INTRODUCTION 


Homopolar  machines  produce  torque  by  the  direct  interaction  of  current,  in  a 
conductor,  with  a  separately  produced  magnetic  flux  in  the  active  region  of  the  machine . 
The  two  primary  functions  of  the  active  region  are  the  transport  of  the  load  current  and 
the  magnetic  flux  at  right  angles  to  each  other.  In  a  homopolar  machine  with 
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superconductive  excitation,  the  conductor  of  choice  to  conduct  the  load  current  through 
the  active  region  of  the  machine  is  copper.  Copper  has  good  electrical  properties  but  has 
a  magnetic  permeability  the  same  as  air.  The  low  permeability  is  not  a  problem  since 
most  of  the  flux  path  is  already  air  and  the  superconductive  magnet  provides  the  required 
ampere  turns  in  a  compact  coil  with  low  power  requirements. 

However,  in  a  homopolar  machine  with  normal  excitation,  a  complete  ferrous 
path  must  be  provided  because  of  the  limited  ampere  turns  available  from  a  normal 
magnet.  In  this  case,  transporting  the  current  through  the  active  region  of  the  machine 
becomes  more  difficult  and  involves  more  tradeoffs.  The  use  of  continuous  copper 
sections  provides  an  excellent  current  path  but  introduces  a  large  effective  air  gap  in  the 
magnetic  circuit.  This  results  in  significant  increases  in  the  size  and  power  consumption 
of  the  field  coils  with  direct  negative  impacts  on  the  efficiency  and  size  of  the  machine. 
One  alternative  is  the  use  of  a  slotted  iron  rotor  with  copper  bars  in  the  slots  to  carry 
current  and  the  ferrous  material  remaining  between  the  slots  to  carry  the  magnetic  flux. 
The  primary  disadvantage  of  this  approach  is  that  the  active  region  must  be  enlarged, 
since  only  a  portion  of  it  is  available  for  each  of  its  two  functions.  A  second  alternative  is 
the  use  of  iron  conductors  to  carry  both  the  load  current  and  the  magnetic  flux.  The 
advantage  of  this  approach  is  that  the  entire  volume  of  the  active  region  is  available  for 
both  flux  and  current  transport  since  the  conductors  are  ferrous.  Because  of  this 
combination  of  functions  the  active  region  is  smaller,  thereby  reducing  the  overall 
machine  size.  This  is  the  approach  that  was  selected  for  the  demonstration  contra  rotating 
motor  being  developed  as  part  of  this  task. 

There  are,  however,  interactions  between  the  load  current  and  the  transverse 
magnetic  flux  that  must  be  analyzed  and  accounted  for  in  the  design  of  the  machine.  The 
load  currents  in  the  iron  bars  generate  an  internal  magnetic  field  intensity  that  must  be 
added  vectorialy  to  the  torque  producing  radial  field  from  the  central  field  coil.  This 
effect  is  shown  in  Fig.  1  where  0a  is  flux  from  the  load  current  passing  axially  through 


2 


the  iron  bars,  and  <})r  is  the  radial  flux  from  the  field  coil.  In  superconducting  machines 
which  use  copper  bars,  <|>a  is  very  small  because  the  permeability  of  copper  is  very  low. 
However,  in  a  normally  conducting  machine  with  iron  bars,  0a  is  significant  and  the 
resulting  combination  of  the  two  fluxes  (shown  as  part  c  of  Fig.  1)  produces  a  higher  total 
magnetic  field  in  the  bar.  All  magnetic  materials  are  non-linear  and  most  have  a  lower 
permeability  at  higher  magnetic  field  intensities.  Thus  the  working  flux  in  the  machine 
sees  a  lower  effective  permeability  in  the  iron  bars  due  to  presence  of  the  axial  current, 
which  must  be  accounted  for  by  increasing  the  amp-tums  of  the  field  coil. 


Fig.  1.  Components  of  Load  Current  Flux  Interaction. 

a)  Circular  flux  from  armature  current  in  upper  and  lower  iron  bars. 

b)  Radial  flux  from  the  field  winding  with  no  armature  current. 

c)  Interaction  of  field  and  armature  flux  when  both  field  and  bar  currents  are  present. 
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A  second  effect  which  may  need  to  be  accounted  for  is  the  increase  in  resistivity  of 
the  iron  conductors  due  to  magneto  resistive  effects.  That  is,  in  the  presence  of  high 
magnetic  fields,  some  irons  tend  to  show  an  increase  in  resistivity.  In  most  machines  this  is 
not  a  problem  because  iron  is  not  normally  used  to  carry  current.  However  in  the  normal 
conducting  homopolar  motor,  an  increase  in  resistivity  of  the  current  carrying  iron  bars  will 
increase  heat  production  and  lower  the  efficiency  of  the  machine.  This  effect  is  probably 
small  and  was  not  modeled.  Rather,  iron  bars  of  the  same  type  as  will  be  used  in  the  demo 
motor  were  experimentally  tested  and  the  magnetic  field  was  shown  to  have  little  effect. 

APPROACH 

The  approach  taken  for  the  flux  -  load  current  interaction  portion  of  this  work  was 
to  develop  analytical  methods  that  could  be  applied  to  homopolar  machine  design  in 
general  rather  than  an  approach  that  addressed  the  problem  only  for  the  particular 
geometry  of  the  motor  under  development.  The  basic  method  we  selected  to  account  for 
the  flux  -  load  current  interaction  was  a  combination  of  two  2  dimensional  finite  element 
models.  The  flux-load  interaction's  influence  on  the  flux  level  in  the  machines  is  actually 
a  3  dimensional  effect.  However,  3  dimensional  models  are  large  and  cumbersome,  and 
each  model  would  be  tied  to  the  specific  geometry  of  the  bars.  Instead,  this  technique 
consisted  of  using  a  2  dimensional  finite  element  program  to  model  the  flux-load 
interaction's  effect  on  the  permeability  of  the  iron  bars.  These  modified  BH  curves  were 
then  used  to  represent  the  bar  iron  in  another  2  dimensional  model  which  calculated  the 
total  working  flux  and  torque  in  the  machine. 

The  capability  of  finite  element  techniques  to  accurately  model  the  flux  load 
interaction  was  first  verified  experimentally  with  an  apparatus  that  could  be  modeled  with 
our  existing  2D  finite  element  software  (PE2D  by  Vector  Fields).  The  geometry  selected 
consisted  of  a  nominal  1  in.  x  3  in.  steel  bar  placed  between  the  poles  of  a  water  cooled 
electromagnet  which  established  a  transverse  magnetic  induction.  At  the  same  time  a 
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current  was  established  in  the  longitudinal  direction  in  the  bar.  The  current  in  the  magnet 
was  linearly  ramped  up  while  the  magnet  current  and  the  average  transverse  magnet 
induction  were  recorded  on  an  X  Y  plotter.  This  was  done  for  bar  currents  of  0,  2,000, 
4,000, 6,000,  8,000,  and  10,000  amps  and  compared  to  the  results  of  a  finite  element 
model  of  the  apparatus. 

Once  the  finite  element  method  was  verified  using  a  single  bar,  the  complex  flux 
interaction  of  multiple  bars  side  by  side,  shown  in  Fig.  2,  was  modeled  using  finite 
element  techniques.  This  finite  element  model  is  more  complicated  because  the  boundary 
conditions  around  any  given  bar  (symmetry  plus  a  constant)  do  not  fit  the  typical 
boundary  conditions  allowed  for  in  most  finite  element  packages.  An  in-house  finite 
element  program  was  written  to  allow  for  this  boundary  condition.  The  in-house  program 
was  verified  by  running  iron  bar  models  using  regular  symmetry  conditions  and 
comparing  the  results  with  the  PE2D  finite  element  package  used  to  model  the 
experimental  test  setup.  After  verification,  the  in-house  program  was  then  run  using  the 
symmetry  plus  a  constant  boundary  condition  for  various  states  of  radial  flux  and  axial 
current  to  create  effective  BH  characteristics.  These  BH  curves  were  then  used  for  the 
iron  bar  material  in  the  2D  finite  element  model  of  the  entire  machine. 


Fig.  2.  Several  Bars  Side  by  Side  to  Show  Complex  Flux  Interaction 
and  Boundary  Conch tion. 
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The  approach  for  the  magneto  resistive  effect  was  experimental  in  nature.  The 
resistance  of  the  steel  bar  was  measured  under  realistic  magnetic  induction  and  current 
conditions  to  determine  if  the  resistance  changed  in  response  to  the  transverse  magnetic 
induction. 

EXPERIMENTAL  DETERMINATION  OF  FLUX-CURRENT  INTERACTION 
TEST  SETUP 

The  basic  test  setup  consisted  of  a  ferrous  bar  with  a  cross  section  of  .972  in.  x 
2.883  in.  This  bar  was  fitted  with  a  trapezoidal  adapter  bar  on  each  edge  to  provide  a 
gradual  transition  to  the  magnet  pole  face.  Four  G-10  (glass-epoxy)  dowels  maintained 
alignment  between  the  test  bar  and  each  adapter.  A  0.045  in.  air  gap  was  maintained 
between  the  test  bar  and  each  adapter  by  plastic  shim  washers  on  the  dowel  pins.  This 
gap  allowed  insertion  of  a  Hall  effect  probe  to  measure  the  transverse  magnetic  induction. 
Figure  3  shows  the  design  of  the  test  bar  and  adapters.  The  test  bar  was  fitted  with 
thermocouples,  voltage  taps,  and  field  pickup  coils.  This  assembly  was  inserted  between 
the  poles  of  a  15  in.  water  cooled  electromagnet  and  connected  to  a  pair  of  5,000  amp 
power  supplies  with  water  cooled  power  cables.  The  instrumented  test  bar  is  shown  in 
place  between  the  poles  of  the  magnet  in  Fig.  4.  The  electromagnet  was  connected  to  a 
programmable  power  supply  controlled  by  an  arbitrary  wave  form  generator.  The 
instrumentation  set  up  is  shown  in  Fig.  5. 
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Fig.  4.  Instrumented  Test  Bar  Installed  Between  the  Magnet  Poles. 
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Fig.  5.  Instrumentation  Setup 
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The  average  value  of  the  transverse  magnetic  induction  was  obtained  by 
integrating  the  output  voltage  from  the  pickup  coils  with  a  low  drift  integrator.  The 
voltage  induced  in  the  field  pickup  coil  by  the  changing  magnetic  flux  is  given  by 

dC> 

v  =  NdT- 


(i) 


Where  <I>  is  the  transverse  magnetic  flux  and  N  is  the  number  of  turns  in  the  pickup  coil. 
This  voltage  is  applied  to  the  input  of  a  voltage  integrator  whose  output  voltage  is  given 
by 


(2) 


Substituting  Eq  1  into  Eq  2  results  in 

i  r  d<i> 

Vout=RcJNdTdt 

0 

which  when  integrated  gives 

vout  =^4>(t)-4>(o)) 

The  integrator  output  is  initially  set  to  0  which  results  in  an  output  voltage  of 

v 

v°ut-  RC 


(3) 


(4) 


(5) 


where  <t»(t)  is  the  change  in  flux  since  the  integrator  output  was  set  to  zero. 

The  relationship  between  magnetic  induction  and  flux  through  a  surface  S  is 
given  by 


O 


(6) 


With  the  geometry  of  this  setup,  ^  is  normal  to  S,  and  since  we  are  interested  in  the 
average  value  of  I?,  Eq  6  becomes 
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O  =  BavgLW  . 


(7) 


Where  Bavg  is  the  average  normal  magnetic  induction  and  L  and  W  are  the  length  and 
width  of  the  surface  S.  Combining  Eqs  5  and  7  and  solving  for  Bavg  we  have 

d  _  VputRC 
BaVg-  nlw  • 


(8) 


Values  of  R  and  C  were  chosen  such  that 


r  _  Vout 

oavg  -  2 


(9) 


Temperatures  from  the  six  type  T  thermocouples  on  the  bar  were  directly 
monitored  and  recorded  by  a  data  acquisition  system.  The  magnet  current  was  monitored 
with  an  in  line  shunt  and  recorded  on  a  X  Y  plotter.  The  bar  current  was  monitored  with 
in  line  shunts  in  the  power  supplies  and  also  recorded  by  the  data  acquisition  system.  The 
voltage  drop  in  the  center  15  inches  of  the  bar  was  monitored  with  6  sets  of  voltage  taps, 

3  sets  on  each  side.  These  were  measured  and  averaged  by  the  data  acquisition  system. 

During  initial  checkout  and  calibration,  a  problem  with  the  setup  was 
encountered.  If  the  magnet  current  was  ramped  up  very  slowly  the  curve  of  integrator 
output  verses  magnet  current  appeared  proper  and  followed  the  same  path  both  during 
ramp  up  and  ramp  down.  If  the  ramp  rate  was  increased,  the  two  curves  diverged  from 
each  other  resulting  in  an  open  loop.  The  ramp  up  curve  was  below  the  slow  ramp  curve 
and  the  ramp  down  curve  was  above  it,  as  illustrated  in  Fig.  6.  The  faster  the  ramp  rate 
the  more  open  the  curves  became.  We  considered  using  a  slow  ramp  rate  but  encountered 
two  secondary  problems.  The  first  was  drift  of  the  integrator  output  voltage  due  to  input 
offset  voltage.  The  second  was  heating  of  the  uncooled  test  bar,  particularly  at  the  higher 
currents.  The  basic  problem  with  the  faster  ramp  had  to  be  solved  to  insure  accurate  data. 

The  most  likely  cause  appeared  to  be  eddy  currents  in  the  copper  cooling  plates  or 
in  the  iron  yoke  of  the  magnet  assembly.  The  copper  cooling  plates  were  the  most  likely 
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source  of  the  problem  since  they  appear  as  low  resistance  shorted  turns  around  the 
magnet  pole  and  link  100  %  of  the  magnet  flux.  The  effective  ampere  turns  in  the 
magnet 

is  equal  to  the  algebraic  sum  of  the  ampere  turns  in  the  coils  plus  the  ampere  turns  due  to 
eddy  currents.  During  ramp  up  the  eddy  currents  reduce  the  effective  ampere  turns  and 
during  ramp  down  they  increase  the  effective  ampere  turns.  This  explanation  is 
consistent  with  our  observations. 

We  decided  to  use  a  fast  ramp  rate  and  to  compensate  for  the  eddy  currents  in  the 
magnet  assembly.  We  know  that  the  induced  voltages  causing  the  eddy  currents  are 
given  by 


V 


dO 

~  dT  • 


(10) 


The  eddy  currents,  I,  are  equal  to 


(11) 


where  V  is  the  induced  voltage  and  R  is  the  effective  resistance  of  the  eddy  current  path  or 

paths.  If  Eq  10  is  substituted  in  Eq  1 1  we  have 

.  _d$  J_ 

1 "  dt  R  ' 


(12) 
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Fig.  6.  Fast  and  Slow  Ramp  Rate  Curves  of  Integrator  Output  Versus  Magnet  Current 
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In  this  equation  dO/dt  is  the  only  parameter  that  changes  with  ramp  rate,  R  is  a  function 
of  the  geometry  and  material  properties  but  would  be  difficult  to  determine  accurately. 
Since  dd>/dt  can  be  directly  monitored  with  a  coil  that  links  the  flux  we  can 
experimentally  compensate  for  the  eddy  currents  without  actually  determining  their 
magnitudes.  A  three  (3)  turn  coil  was  placed  around  the  electromagnet  pole  and  its 
output  connected  to  a  10  turn  potentiometer  used  as  a  voltage  divider.  The  output  from 
the  voltage  divider  was  connected  in  series  with  the  output  from  the  magnet  current  shunt 
and  the  series  combination  was  connected  to  the  X  Y  plotter.  Thus  a  voltage  term 
proportional  to  dd>/dt  was  subtracted  from  the  current  signal.  By  trial  and  error  a 
potentiometer  setting  was  found  that  resulted  in  good  closed  curves  of  flux  verses  magnet 
current.  These  curves  did  not  change  with  ramp  rate  and  they  matched  the  curves 
obtained  at  a  slow  ramp  rate.  This  setting  was  used  during  all  subsequent  measurements. 
This  corrected  magnet  current  is  referred  to  as  the  effective  magnet  current. 

Magnetic  measurements  are  normally  made  on  demagnetized  samples  to  eliminate 
the  effect  of  residual  magnetization  on  the  results.  In  this  setup  it  was  not  practical  to 
demagnetize  the  sample;  therefore  a  compensating  method  was  developed.  This 
consisted  of  making  a  first  set  of  measurements  with  the  magnetic  induction  in  the  same 
direction  as  the  residual  magnetization  in  the  test  bar  and  then  reversing  the  input  to  the 
magnet  and  repeating  the  measurements  with  the  magnetic  induction  in  the  opposite 
direction.  These  two  sets  of  measurements  were  averaged,  as  shown  in  Fig.  7,  giving  the 
same  results  as  if  the  sample  had  been  demagnetized  prior  to  each  set  of  measurements. 

PROCEDURE 

The  axial  current  in  the  test  bar  was  first  set  to  the  desired  value  with  the  power 
supplies  in  the  current  regulating  mode  to  prevent  drift.  The  integrator  output  was  set  to 
zero  and  the  current  in  the  magnet  was  then  ramped  up  to  a  predetermined  level  and 
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returned  to  zero,  at  a  constant  rate,  by  means  of  the  wave  form  generator  programmed 
with  a  triangular  wave.  The  outputs  from  the  integrator  and  the  magnet  current  shunt 
were  recorded  on  a  XY  potter  as  the  current  was  ramped  up  and  then  down.  The  magnet 
leads  were  then  reversed  and  the  above  procedures  repeated  exactly  with  the  results 
recorded  on  the  same  graph  in  the  XY  plotter. 

The  test  bar  was  allowed  to  cool  and  the  measurements  were  repeated  for  each  of 
the  bar  currents.  The  curves  were  then  digitized  to  allow  averaging  of  the  two  curves  for 
each  test  condition.  The  vertical  axis  was  scaled  to  indicate  Bavg»  the  average  transverse 
magnetic  induction. 
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Fig.  7.  Typical  Curves  of  Integrator  Output  Versus  Magnet  Current 
and  Average  of  the  Two  Curves 
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RESULTS 


The  results  for  each  of  the  bar  test  currents  are  shown  in  Figs.  7  through  13. 

The  experimental  curves  shown  are  the  averages  of  the  two  sets  of  measurements  as 
previously  discussed.  The  discrete  data  points  shown  are  the  results  from  the  finite 
element  model  of  the  experiment.  As  expected,  the  axial  current  in  the  bar  had  a 
significant  effect  on  the  magnet  current  required  to  produce  a  given  average  transverse 
magnetic  induction.  It  is  interesting  to  note  the  trend  for  all  the  curves  to  converge  at 
high  values  of  magnet  current,  as  shown  in  Fig.  14,  which  corresponds  to  high  values  of 
magnetic  field  intensity.  At  high  values  of  magnetic  field  intensity  the  steel  is  at  or  near 
magnetic  saturation  and  the  addition  of  the  bar  current  has  little  effect.  It  should  also  be 
noted  that  despite  the  similarity,  these  curves  are  not  the  same  as  normal  magnetization 
curves  for  steel.  The  two  air  gaps  and  the  iron  in  the  electromagnet  which  are  included  in 
the  path  have  a  significant  impact  on  the  results.  In  addition,  the  effect  of  the  load  current 
on  the  effective  magnetic  properties  is  dependent  on  the  geometry  of  the  bar.  Thus  the 
results  do  not  represent  material  properties  as  such,  but  are  intended  to  verify  the  finite 
element  modeling  approach  as  discussed  in  the  next  section. 
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vg.  Transverse  Magnetic  Induction  (Tesla) 


Fig.  9.  Magnetic  Induction  Venus  Magnet  Current 
at  2,000  Amps  Bar  Current 
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vg.  Transverse  Magne 
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Fig.  10.  Magnetic  Induction  Versus  Magnet  Current 
at  4,000  Amps  Bar  Current 
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Fig.  11.  Magnetic  Induction  Versus  Magnet  Current 
at  6,000  Amps  Bar  Current 
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vg.  Transverse  Magnetic  Induction  (Tesla 


0  2  4  6  8  10  12  14 


Effective  Magnet  Current  (amps) 


Fig.  13.  Magnetic  Induction  Versus  Magnet  Current 
at  10,000  Amps  Bar  Current 
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Fig.  14.  Magnetic  Induction  Versus  Magnet  Currv  t 
at  0, 2, 4, 6, 8,  and  10  Kamps  Bar  Current 
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VERIFYING  FINITE  ELEMENT  METHOD 

As  is  evident  from  the  curves  shown  in  Figs.  8  through  13,  the  finite  element 
model  results  correlate  very  well  with  the  experimental  data.  For  all  cases  of  bar  currents 
less  than  6  K  amps  or  field  currents  less  than  6  amps,  the  variations  between  the 
measured  and  calculated  values  are  within  experimental  error  and  are  not  significant.  The 
error  in  the  cases  where  the  bar  currents  were  more  than  6  K  amps  and  field  currents  were 
more  than  6  amps  is  probably  a  result  of  the  fact  that  the  iron  BH  curve  used  in  the  finite 
element  models  was  not  an  exact  match  for  the  iron  actually  used  in  the  tests.  Overall, 
this  experiment  verifies  that  the  finite  element  approach  that  we  developed  accurately 
models  the  experimental  setup  and  can  be  used  for  the  model  motor  under  development. 
This  finite  element  approach  is  versatile  and  can  be  used  for  other  machines  that  use 
ferrous  conductors  in  the  active  region. 

EXPERIMENTAL  DETERMINATION  OF  MAGNETO  RESISTIVE  EFFECTS 
TEST  SETUP 

The  basic  test  setup  was  the  same  as  previously  described  with  the  exception  that 
some  of  the  measurements  were  done  in  a  different  manner.  The  magnetic  induction  wa„ 
directly  measured  by  means  of  a  flat  Hall  effect  probe  inserted  in  the  gap  between  the  test 
bar  and  the  trapezoidal  adapter,  and  all  other  parameters  were  recorded  by  a  data 
acquisition  system. 

PROCEDURE 

With  the  test  bar  in  position  between  the  poles  of  the  electromagnet,  a 
longitudinal  current  of  2,000  amps  was  established  in  the  bar.  The  voltage  drop  from  the 
voltage  taps  and  the  temperatures  were  recorded  by  the  data  acquisition  system  as  the  bar 
was  allowed  to  heat  to  approximately  170  °F.  The  bar  current  was  reduced  to  zero  and 
the  bar  was  cooled  to  room  temperature.  A  transverse  magnetic  induction  of  0.5  Tesla 
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was  then  established  by  the  electromagnet.  The  longitudin  il  current  of  2,000  amps  was 
reestablished  and  the  voltage  and  temperature  measurements  were  again  recorded  as  the 
bar  heated.  This  procedure  was  repeated  for  magnetic  inductions  of  1.0, 1.5,  and  2.0 
Tesla  at  both  2,000  and  4,000  ampere  longitudinal  currents.  Voltage  drops  were 
extracted  from  each  of  the  data  sets  for  temperatures  of  100,  120,  140,  and  160  °F.  The 
six  voltage  drops  for  each  test  condition  were  averaged  and  the  bar  resistance  calculated. 

RESULTS 

The  test  results  for  2,000  amps  are  shown  in  Fig.  15  and  for  4,000  amps  in 
Fig.  16.  As  is  evident  from  the  curves,  there  is  no  significant  variation  in  resistance  of  the 
steel  bar  as  a  result  of  the  transverse  magnetic  induction.  The  maximum  variation  between 
the  2,000  and  4,000  amp  data  is  0.7  %  and  the  maximum  change  with  magnetic  induction  is 
0.4  %.  Both  of  these  are  well  within  experimental  error  for  this  setup.  These  variations  are 
not  significant  in  the  design  of  a  machine  and  do  not  need  to  be  accounted  for. 


Fig.  15.  Resistance  Versus  Magnetic  Induction  at  2,000  Amperes. 
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Fig.  16.  Resistance  Versus  Magnetic  Induction  at  4,000  Amperes. 


FINITE  ELEMENT  ANALYSIS  OF  FLUX  LOAD  INTERACTION 
PROCEDURE 

The  complex  flux  interaction  of  multiple  bars  side  by  side  (shown  Fig.  2)  was 
modeled  using  finite  element  techniques.  This  finite  element  model  is  more  complicated 
because  the  boundary  conditions  around  any  given  bar  do  not  fit  the  typical  boundary 
conditions  allowed  for  in  most  finite  element  packages.  A  close  examination  of  Fig.  2 
shows  that  there  is  a  symmetry  condition  between  each  bar.  That  is,  the  potential 
functiw.  used  to  describe  one  bar  can  describe  any  bar  because  the  flux  patterns  are  the 
same.  The  potential  function  for  each  bar  is,  however,  not  numerically  equal.  The 
uniform  radial  field  implies  that  there  must  be  a  potential  gradient  circumferentially 
through  the  airgap  so  each  bar  will  be  at  a  higher  potential  level  than  the  bar  next  to  it. 
The  complete  radial  field  is  three  dimensional,  therefore,  it  is  impossible  to  model  all  the 
bars  in  the  machine  at  once  in  only  two  dimensions  because  the  radial  field  would  apear 
from  a  point  (not  mathematically  possible).  However,  a  single  upper  and  lower  pair  can 
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be  modeled  in  two  dimensions.  It  can  be  seen  by  the  symmetry  of  Fig.  2  that  an 
appropriate  boundary  condition  for  any  pair  of  bars  (upper  and  lower)  is  that  for  any 
given  height,  the  magnetic  potential  at  the  center  of  the  air  space  to  the  left  of  the  bars  is 

constrained  to  be  equal  to  the  potential  at  the  center  of  the  gap  to  the  right  of  the  bars  plus 
a  constant  to  allow  for  the  radial  flux  0r  This  symmetry  plus  a  constant  boundary 

condition  was  not  available  in  the  finite  element  package  used  to  model  the  experimental 
tests,  so  an  in-house  finite  element  program  was  written  to  allow  for  this  boundary 
condition.  The  theory  used  in  the  in-house  program  is  described  in  Appendix  A,  and  a 
listing  of  the  program  appears  in  Appendix  B. 

Figure  17  is  a  scale  drawing  of  the  finite  element  mesh  used  to  model  the  flux-load 
current  interaction.  It  consists  of  two  bars  with  a  solid  iron  shield  above  and  solid  iron  core 
below.  The  dimensions  of  the  bar  model  are  shown  in  Fig.  18.  Each  triangular  element  of 
the  model  and  the  node  comers  of  all  elements  are  shown.  The  currents  in  the  two  bars  are 
going  in  opposite  directions.  The  symmetry  plus  a  constant  boundary  condition  is  applied 
at  the  center  of  the  space  between  the  bars.  The  output  of  the  finite  element  solution  is  the 
magnetic  potential  at  each  node,  from  which  B  and  H  data  can  be  calculated  for  the  given 
radial  flux  and  current  condition. 
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Vn=V(n.1760)- AM< 


Fig.  17.  Finite  Element  Model  Showing  the  Elements  Used  and  Boundary  Conditions 
for  the  Case  Where  the  Armature  Current  is  35  kA.  The  Comer  of  Each 
Triangular  Element  Represents  a  Node,  Numbered  1  to  1840  Starting 
From  the  Lower  left  Hand  Comer. 
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Fig.  18.  Dimensions  of  Finite  Element  Bar  Model. 
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CALCULATIONS 

One  parameter  that  can  affect  the  accuracy  of  the  solution  the  mesh  density.  The  smaller 
each  individual  element,  the  more  accurate  the  solution,  until  numerical  roundoff  becomes  an 
issue.  For  the  computer  for  which  the  program  was  written,  a  model  of  1 840  nodes  and  1738 
elements  (700  nodes  per  meter)  was  as  large  as  could  be  comfortably  run  without  memory 
problems.  Figure  19  shows  some  JVfrdl  integration  results  for  various  node  densities,  for  an  iron 
bar  run  with  B-radial  =  .76  Tesla  and  I-bar  =  35000  Amp.  By  extrapolating  Fig.  19,  it  is 
estimated  that  a  node  density  of  700  nodes/meter  will  produce  about  a  7%  mesh  error,  which  is 
reasonable  considering  the  errors  inherent  in  measuring  consistent  BH  properties. 

Many  iron  bar  finite  element  models  were  run  using  radial  flux  densities  from  0.1  to  8 
Tesla,  and  armature  currents  from  0  to  100,000  Amps.  Some  sample  flux  plots  are  shown  in 
Fig.  20  for  B-radial  =  .76  Tesla  and  various  current  densities  in  the  bars.  The  flux  plots  in 
Fig.  20  show  lines  of  constant  potential,  the  density  of  lines  being  equivalent  to  flux  density. 
Note  the  familiar  fingerprint  pattern  similar  to  Fig.  2.  Also  note  that  the  flux  densities  in  the 
bars  increases  significantly  even  though  the  flux  densities  in  the  core  and  shield  stay  constant. 
The  next  step  was  to  calculate  the  equivalent  BH  characteristics  for  the  iron  bars  which 
involves  determining  the  actual  H  that  exists  in  the  iron  for  any  given  B-radial  seen  by  the 
homopolar  motor  for  a  particular  axial  current  case.  The  equivalent  BH  characteristics  were 
needed  for  an  already  existing  finite  element  model  of  the  overall  homopolar  machine.  The 
overall  model  included  the  airgaps  above  and  below  the  bars,  but  not  the  radial  spacers  (which 
would  require  3  dimensions).  Therefore  the  airgaps  are  not  modeled  correctly  because  the 
iron  bars  concentrate  the  flux  underneath  and  above  them,  so  that  Jh»(11  through  the  airgaps  is 

larger  than  it  would  be  if  the  flux  were  evenly  distributed  over  the  airgap.  This  was  accounted 
for  by  modifing  the  values  of  Hfe  obtained  from  the  flux  bar  model. 
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Bar 

Current 

Density 


B.  Rad’dl 


Fig.  20.  Flux  Plots  in  Ferrous  Conductors  with  External  Magnetic  Fields. 

Contours  Represent  Lines  of  Constant  Potential. 

Contour  Density  is  Equivalent  to  Flux  Density. 
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Average  values  for  magnetic  field  intensity  H  can  be  obtained  from  the  flux  bar 
model  by  integrating  H»dl  along  a  line  sement  and  then  dividing  by  the  length  of  the 
segment.  Using  this  method  the  equivalent  magnetic  intensity  Hfe  for  the  iron  bars  was 

calculated  from 


H  •  dl 


a 


•irgaps 


B 


radial 

^0 


dl 


'fe 


(13) 


where 

r 

H  •  dl  is  integrated  from  the  top  of  the  rotor  core  to  the  bottom  of  the  stator 

tout 


shield 

|T®radial 
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•irgaps 


*0 


J 


dl 


is  integrated  through  the  bottom,  middle,  and  upper  airgaps 


and  Lfe  is  the  total  height  of  the  iron  bar  segments. 


Note  that 


Jfl*dl  is 


the  total  Amp-turns  needed  to  cross  the  iron  bars  and  airgaps  under 


total 

actual  conditions,  and 


•irgaps 


^radial  ^ 

V  "0  , 


B, 


radial 


dl  =  ^airgaps  is  the  Amp-turns  needed  for  an 


evenly  distributed  radial  flux  to  cross  the  airgaps  only  which  is  already  included  in  the 
overall  model.  A  plot  of  Hfc  verses  B-radial  for  all  cases  is  shown  in  Fig  19  and  Table  1. 

This  is  the  data  that  was  used  the  homopolar  motor  overall  model. 
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RESULTS 


The  plots  in  Fig.  21  and  the  data  in  Table  1  represent  the  final  product  for  this 
segment  of  the  project.  These  curves  show  that  the  effective  permeability  of  the  iron  in 
the  bars  is  significantly  affected  by  the  presence  of  axial  current.  These  curves  were  used 
to  represent  the  bars  in  a  finite  element  model  of  the  normal  conducting  homopolar  demo 
motor.  The  motor  analysis  will  be  described  in  a  separate  report.  However,  in  the  case  of 

the  demo  motor,  the  effect  of  the  flux  concentration  due  to  axial  currents  turned  out  to  be 
less  than  expected.  In  the  nominal  full  power  case  of  1^  =  35000  A  and  BradiaJ  =  .76 

Tesla,  the  modified  BH  curves  only  reduced  the  overall  working  flux  and  torque  of  the 
machine  by  4%.  This  is  probably  because  the  iron  bars  are  only  a  short  segment  of  the 
path  of  the  working  flux,  and  much  of  the  rest  of  that  path  includes  iron  near  saturation. 


35 


Fig.  21 

Effective  BH  Characteristics  for  Bar  Rotor  with  Indicated  Armature  Current 


Note  to  page  21  -  Effective  BH  characteristics  for  iron  bars  with  axial  current  include  flux 
concentration  due  to  radial  spacers  between  the  bars  (Hfe).  The  iron  bars  and  radial 
spacers  are  treated  as  a  single  material,  therefore  this  data  is  valid  only  for  the  geometry 
modeled.  These  curves  were  used  to  represent  the  iron  bars  in  a  finite  element  model  of 
the  entire  homopolar  machine. 
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Table  1.  Data  for  BH  Curves  in  Fig.  21 


Bar  Current 

B  Radial 

H  Iron+Spacers 

(A) 

CD 

(At/m) 

OK 

0.100 

OK 

0.400 

908.0 

OK 

0.760 

1730.0 

OK 

0.875 

1990.0 

OK 

1.000 

2300.0 

OK 

1.200 

2850.0 

OK 

1.400 

3830.0 

OK 

1.600 

9020.0 

OK 

1.800 

21500.0 

OK 

2.000 

55900.0 

OK 

8.000 

4730000.0 

35K 

0.100 

739.0 

35K 

0.400 

3000.0 

35K 

0.760 

5640.0 

35K 

0.875 

6640.0 

35K 

1.000 

7610.0 

35K 

1.200 

9640.0 

35K 

1.400 

11800.0 

35K 

1.600 

14400.0 

35K 

1.800 

24000.0 

35K 

2.000 

60900.0 

35K 

8.000 

4730000.0 

67.5K 

67.5K 

0.100 

0.400 

67.5K 

0.760 

8970.0 

67.5K 

0.875 

10300.0 

67.5K 

1.000 

11900.0 

67.5K 

1.200 

15000.0 

67.5K 

1.400 

18300.0 

67.5K 

1.600 

22600.0 

67.5K 

1.800 

30700.0 

67.5K 

2.000 

72700.0 

67. 5  K 

8.000 

4730000.0 

100K 

0.100 

1780.0 

100K 

0.400 

6600.0 

100K 

0.760 

12300.0 

100K 

0.875 

13800.0 

100K 

1.000 

16200.0 

100K 

1.200 

20400.0 

100K 

1.400 

25000.0 

100K 

1.600 

30900.0 

100K 

1.800 

41700.0 

100K 

2.000 

88200.0 

100K 

8.000 

4730000.0 
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CONCLUSIONS 


An  accurate,  experimentally  verified  model  of  the  flux-load  current  interaction 
has  been  devised.  It  shows  that  the  permeability  of  iron  bars  changes  significantly  in  the 
presence  of  axial  currents,  and  that  this  effect  must  be  taken  into  account  for  an  accurate 
calculation  of  machine  flux  and  torque.  The  model  is  flexible,  and  can  be  used  for  other 
bar  geometries.  Although  developed  for  a  normal  conducting  homopolar  motor,  the  in- 
house  finite  element  program  can  be  useful  for  other  configurations  that  require  the 
special  "symmetry  plus  a  constant"  boundary  condition,  such  as  modeling  a  single  slot  of 
an  AC  machine. 
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APPENDIX  A 


FINITE  ELEMENT  THEORY 

Finding  the  magnetic  field  produced  by  electric  currents  is  straightforward  as  long 
as  the  current  and  magnetic  field  exist  entirely  in  a  uniform  linear  medium  such  as  air. 

But  once  material  boundaries  and  non-linear  materials  are  introduced,  magnetic  field 
problems  often  become  too  complex  to  be  solved  by  direct  analytical  techniques.  Such 
cases  often  use  the  finite  element  method,  in  which  a  complex  combination  of  materials 
and  geometry  is  modeled  by  many  small  elements,  each  consisting  of  a  uniform  material 
of  constant  permeability.  The  large  problem  is  then  solved  by  finding  the  boundary 
conditions  that  will  solve  each  of  the  small  element  problems  simultaneously.  Any  non¬ 
linear  effects  where  the  permeability  of  a  material  is  dependant  on  the  solution  are  then 
resolved  by  iteration. 

The  first  step  is  to  find  an  equation  relating  the  magnetic  potential  boundary 
conditions  around  a  small  element  of  uniform  material  to  the  current  passing  through  the 
element.  In  two  dimensions,  the  magnetic  potential  ^  can  be  found  from  Poisson’s 
general  field  equation 


'-h  ♦  *  a 

p  o x  m  oy 


0 


(A  1.1) 


were  Q  is  current  density,  p.  is  magnetic  permeability,  and  x  and  y  are  Cartesian 
coordinates.  Unlike  most  potential  problems,  Q  and  \J/  are  actually  vector  quantities. 
However,  in  two-dimensional  analysis,  the  direction  of  Q  and  t/r  are  always  in  the  z 
direction  (perpendicular  to  x  and  y )  so  that  numerically  they  can  be  treated  a  scalar 
quantities. 


To  apply  this  equation,  assume  we  have  a  triangular  element  of  a  uniform  linear 
material  with  three  nodes  i,  j,  and  k.  Further,  we  assume  that  t/r  varies  linearly  over  the 
element,  which  is  a  reasonable  assumption  as  long  as  the  element  is  small  compared  to 
the  entire  problem.  The  node  potential  equations  for  a  linear  triangular  element  that 

satisfy  Poisson’s  equation  are  found  in  Segerlind1.  The  form  of  the  potential  function  is 

[K]{lf}  -  (F)  =  {0}  (A  1.2) 

where  K  is  the  “stiffness”  matrix,  iff  is  the  scalar  potential  at  each  node,  and  F  is  the 
“forcing  function.”  The  forcing  function  is  calculated  as 

{ F }  =  -Op '  1  ’ 

J  1 1 J 


1  Larry  J.  Segerlind,  Applied  Finite  Element  Analysis, {2nd  Edition)  Wiley:  New  York,  1984  pp  56-58, 9 1  -93, 
115-120  and  165-171.  TA  347  F5  S43  1984 
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The  stiffness  matrix  is  calculated  as 


[K] 


1 

4/j.A 


bf 


bfij 


bih 


bibj 

bibk 

bj 

bjbk 

bjbk 

1 

4  fiA 


CiCj 


Wk 


CiCj 

cick 

'1 

cjck 

cjck 

4 

where  A  is  the  area  of  the  triangular  element  and  constants  b  and  c  are  calculated  from  the  x 
and  y  coordinates  for  nodes  i,  j,  and  k  as  follows: 

bi  =  y}  -  yk  Ci  =  xk  -  x} 

bj  =  yk  -  y,-  cj  =  Xi  -  xk 

bk  =  yt  -  yj  ck  =  *j  *  *i 

Note  that  the  area  of  a  triangle  is 

A  .  ^(x,yj*xjyk*xkyl-xiyt-xJyrx ky) 


The  area  is  positive  if  the  order  of  the  nodes  is  chosen  so  that  the  surface  vector  is  in  the 
positive  z  direction.  That  is,  if  side  1  x  side  2  and  side  2  x  side  3  point  in  the  z 
direction. 

Assuming  the  permeability  /x  is  constant  for  any  given  iteration,  Eq  A  1.2  is  a  set  of  three 
linear  algebraic  expressions  per  element  relating  the  potential  at  the  three  nodes  to  the 
current  in  the  element  and  the  element  shape.  When  elements  are  side  by  side  there  will 
be  more  than  one  equation  per  node.  But  since  these  are  linear  equations,  the  various 
equations  for  each  node  can  be  added  together  so  that  there  will  be  as  many  equations  as 
nodes.  The  potential  of  some  nodes  may  be  set  to  a  constant  or  specified  by  a  given 
function  of  other  nodes  to  allow  for  special  boundary  conditions.  For  example,  in  this 
project,  the  potential  of  the  nodes  at  the  far  right  side  of  the  model  were  set  equal  to  a 
constant  plus  the  potential  of  the  corresponding  nodes  at  the  far  left  side  of  the  model. 

The  potentials  that  are  not  set  to  define  boundary  conditions  are  found  by 

(?)  *  [K] "'(  F) 


were  [  K]  is  found  numerically.  Since  [  K]  is  usually  a  large  matrix,  finding  its  inverse 
involves  much  computation,  and  many  sophisticated  techniques  have  been  derived  for 
efficient  solution.  However,  for  the  purposes  of  this  project,  a  basic  iteration  technique 
was  used  since  computation  speed  was  not  critical  and  it  is  the  easiest  to  program. 

Once  the  potential  at  each  node  is  calculated,  the  flux  density  can  be  calculated  from  the 
derivative  of  the  potential 


B, 


dy 

dy 


and  By 


d\ff 

dx 
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Numerically  this  can  done  by  writing  an  equation  for  the  linear  change  of  iff  over  the 
element  where  the  end  potentials  are  known.  That  is 


iff  =  +  <*2  x  +  03  y 


where  by  definition 


The  boundary  conditions  are 


B . 


dy  ~  a 3 


5, 


dx  "  °2 


%  =  «1  +  °2Xi  *  “s?, 
Vj  =  aj  ♦  c^xj  +  o^y j 

+  <hxk  + 


Solving  for  av  a 2,  and  Oj  produces 


%(yj*y *)  *  %(yk*yi)  *  vfk<yi*yj) 

“2  s  By  =  2A 


%(xk*Xj)  *  %(Xj*xk)  *  %(Xj  +  Xj) 

“3  =  Bx  =  2A 

and  By  are  then  used  to  find  fi  for  each  non-linear  element  using  a  lookup  table  with 
the  BH  characteristics  for  that  material.  The  whole  process  is  then  repeated  until  /x  and  iff 
no  longer  change. 
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APPENDIX  B 


This  appendix  is  a  listing  of  the  finite  element  program  developed  in-house  at  the 
Annapolis  Detachment,  Carderock  Division,  Naval  Surface  Warfare  Center  to  calculate 
flux-current  interactions  with  symmetry  plus  a  constant  boundary  conditions.  The 
program  was  written  in  Microsoft  QuickBasic  to  run  on  a  Macintosh  Hex  computer 
running  system  6.0x.  (The  program  needs  some  minor  revision  to  work  with  the  new 
system  7).  The  program  is  in  five  parts:  1)  user  interface  and  main  program,  2)  model 
creator,  3)  find  node  and  element  connections,  4)  write  equations,  and  5)  solve  equations. 
Each  part  is  listed  separately  with  a  description  of  the  significance  of  each  part. 


t  re  i  dim3dki  ici  i  an  as  a  Kiiei  rcitS 


This  is  the  main  program  from  which  the  other  programs  are  called.  It  includes  a  menu 

driven  interface  which  calls  other  programs  that  create  and  solve  the  model.  Once 

created,  the  model  can  be  examined  both  graphically  and  directly.  Some  examine 

routines  require  information  that  is  not  available  until  the  equation  writer  and  solver  have 

been  called.  A  special  routine  added  to  this  program  is  "Integrate  H*dl",  which  calculates 

J  H*dl  between  specified  points.  This  routine  was  used  to  calculate  the  effective  BH 

characteristics  of  the  iron  with  axial  currents. 

’  Finite  Element  Main  Program  9/15/89 

COMMON  WhoAmI$,FileName$,Printf%,menuO%,menul  % 

WhoAmI$="FE  8/21  shell.s" 

Solve$="Solve  Mesh.s" 

Connect$="Find  Node  Connections.s" 

WriteEqu$="  Write  Equations.s" 

IF  SYSTEM(4)  THEN  *  SYSTEM(4)=Tnje  if  program  is  compiled 
WhoAmI$="" 

SolveS="Solve  Mesh" 

Cdnnect$="Find  Node  Connections" 

WriteEqu$="Write  Equations" 

FOR  adr%=&H91 1  TO  &H910  +  PEEK(&H910)  ’  Looks  up  name  of  Application 
WhoAmI$=WhoAmI$+CHRS(PEEK(adrffc)) 

NEXT  adr% 

END  IF 

PRINT  "**"+WhoAmI$+"**" 

OPTION  BASE  0 

DIM  Nu%(2000),Vc(2000),x(2000),y(2000),v(2000) 

DIM  Ni%(2(X»),Nj%(2000),Nk%(2000),Nl%(2000) 

DIM  mate%(2000),Type%(2000),Curd(2000),ul(2000),u2(2000) 

DIM  Na%(2000),Nb%(2000),Nc%(2000),Nd%(2000) 

DIM  plot(200) 

MENU  1,0,1, "File" 

MENU  1,1, 1, "Open" 
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MENU  13.0,’Save  As" 

MENU  13.1. "Quit" :  CmdKey  1.3."Q" 

MENU  2,0,1. "Create" 

MENU  2,1,1, "Run  External  Program  Create  Mesh" 

MENU  23,0,"-" 

MENU  23.1. "Find  External  Program" 

MENU  2,4,1, "Re-Write  Equations" 

MENU  3,0,1, "Solve" 

MENU  3,1,1, "Solve  Mesh" 

MENU  4,0,l,"Examine" 

MENU  4,1,1,’Examine  Entire  File" 

MENU  4 3,1, "Examine  Node  Positions" 

MENU  43,1, "Examine  Node  Potentials" 

MENU  4,4,l,”Examine  Node  Connections" 

MENU  43,l,"Examine  Element  Compostion" 

MENU  4 ,6,1, "Examine  Element  Properties" 

MENU  4,7,l."Examine  Element  Connections" 

MENU  43,l,"Examine  Solution  Equations" 

MENU  43.0."-" 

MENU  4,10,l,"Print  to  Screen" 

MENU  4,11.1  ."Print  to  File" 

MENU  4,12,l,"Cancel  Printing" 

MENU  5,0,1, "Plots" 

MENU  5. 1,1, "Material  Outline" 

MENU  53.1. "Mesh" 

MENU  53.1,"Regions" 

MENU  5,4,1, "Flux" 

MENU  6,0,1, "Integrate" 

MENU 6,l,l,"H*dl"  :  CmdKey  6,1,"H" 

'MENU  63.1. "print  V70  equations" 

u0=.000001256637# 

MeshFile%=0  ’  Nothing  has  been  read  from  the  file  yet  (even  if  returning) 

CLOSE  #2  '  Always  return  with  screen  as  output  for  examine  menu 

OPEN  "SCRN:"  FOR  OUTPUT  AS  #2  '  May  want  to  eliminate  later 

RetumFn%=l 

PRINT  "***"JFileName$;"***";menuO%^nenul% 

IF  FileNameSo""  THEN 

MulUCommand%=l  :  menu0p%=menu0%  :  menu  1  p%=menu  1  % 

IF  menu0%=2  AND  menu  1  %<4  THEN  MultiCommand%=0 
IF  menu0%=0  THEN  MultiCommand%=0  '  Set  by  Solve  Mesh  prog 
RetumFn%=0 

menu0%=l  :  menul%-l  :  GOTOwhile2 
END  IF 

loop:  menu0p%=0 
menulp%*0 
menu0%=MENU(0) 
mselecL  menul%=MENU(l) 

Multi  Com  mand%=0 
level%=0 

while2:  WHILE  menu0%  OR  MulliCommand% 
scase:  SELECT  CASE  menu0% 

•*•***••*•**  CASE  1  Open  File  **»••**»•*•** 

CASE  1  'File  Menu  selected 
SELECT  CASE  menu  1% 

CASE  1 
’  Selected  Open. 

IF  RetumFn%  THEN  FileNameS=FILES$(l,"MESH") 
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IF  FileNameS=""  THEN 
MuluCommand%=0 
ELSE 

changeCursor  4 
RctumFn%=l 

CLOSE  #1  '  Close  File  if  already  open 

OPEN  "R",#l  FileNameS, 8 
WINDOW  1 ,  FileNameS 

FIELD  #1, 2  AS  aa2S,  2  AS  ab2$,  2  AS  ac2S,  2  AS  ad2S 
FIELD  #1, 2  AS  ba2$,  2  AS  bb2$,  4  AS  bc4$ 

FIELD  #  1 , 4  AS  ca4$,  4  AS  cb4$ 

GET  #1.1 

Nodes%=CVI(aa2S) :  Elements%=CVI(ab2$) 

Nodes  l%=Nodes%+l :  Elements  l%=Elements%+l 
TestV %=CV I(ac2S) :  TestE%=CVI(ad2S) 

PRINT  "Nodes  =";Nodes%;"  Elements  =";Elements% 

PRINT  TestV%,TestE% 

nuvcF%=0 :  inuvc%=l 

xyF%=0 :  ixy%=l  :  xypt&=l+Nodes% 

vF%=0:  iv%=l :  vpt&=l+2*Nodes% 

nijklF%=0  :  inijkl%=  1  :  nijldpt&=l+3*Nodes% 

mcurdF%=0 :  imcurd%=l  :  mcurdpl&=l+3*Nodes%+Elements% 

typeF%=0 :  itype%=l :  typept&=l+3*Nodes%+Elements% 

uF%=0 :  iu%=l :  upc&=l+3*Nodes%+2*Elements% 

nregF%=0:  ireg%=l :  regptA=l+3*Nodes%+3*Elements%  :  nreg%=0 

nabcdF%=0:  inabcd%=l  :  nabcdpt&=l+3*Nodes%+4*Elements% 

equpt&=2+3*Nodes%+5*Elements% 

MeshFile%=l 
MENU  1,2,1, "Save  As" 

END  IF 

CASE  2 

*  Selected  Ssvc  As 

NewMFUeS=FILESS(0."New  MESH  Filename") 

IF  NewMFileSo""  THEN 
changeCursor  4 

copyFile  FileNameS,  NewMFileS 
NAME  NewMFileS  AS  NewMFileS.  "MESH" 

FileNameS=NewMFileS 
WINDOW  1,  FileNameS 
END  IF 

CASE  3 

'  Selected  QuiL  Be  sure  to  save  current  results 
STOP 

case  else 

PRINT  "Error.  Not  Expecting  Menu  1  Item";MENU(l);"to  be  selected." :  INPUT  x :  STOP 
END  SELECT 


•••»••••••*•  CASE  2  Run  Model  Creator  ************* 

CASE  2  '  Create  Menu  selected 
SELECT  CASE  menu  1  % 

CASE  1 

IF  SYSTEM(4)  THEN  mneS=FILESS{l."Pext")  ELSE  mneS=FILESS(  1  ."TEXT") 
IF  mneSo""  THEN  CHAIN  mneS 

CASE  3 

'  Change  File  Type  to  Pext  (external  Program) 

IF  SYSTEM(4)  THEN  mneS=FILESS(l,"APPL")  ELSE  mne$=FHJESS(l,,TEXT") 
IF  mneSo""  THEN 

IF  SYSTEM(4)  THEN  NAME  mneS  AS  mneS,  "Pext" 
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CHAIN  mne$ 

END  IF 

CASE  4 

'  Manual  re-write  equations 
IF  MeshFik%  THEN 
GOSUB  findPnode 

menu0%=2  :  menul%=5  :  GOTO  scase 
ELSE 

PRINT  "No  File  Open." 

END  IF 

CASE  5 

menu0%=2:  menul%=6 
CHAIN  Connects 

CASE  6 

menu0%=2 :  menul%=7 
CHAIN  WriteEquS 

PACE  7 

PRINT  "Completed” 

MultiCommand%=0 
GOTO  Idle 

CASE  ELSE 

PRINT  "Error  Not  Expecting  Menu  2  Item";MENU(l);"to  be  selected." 

INPUT  x:  STOP 
END  SELECT 

•••*•*•**•••  CASE  3  Run  Equation  Solver  ************* 

CASE  3  '  Solve  Menu  selected 

’  IF  TestV%>3  THEN  PRINT  "Convergence  already  completed" :  GOTO  Idle 
changeCursor  4 
IF  MeshFile%=0  THEN 
MultiCommand%=l 
menu0p%=3 
mcnulp%=l 

menuO%=l :  menul%=l :  GOTO  scase 
END  IF 

MultiCommand%=0 

IF  (TestV%  AND  2)=0  THEN  GOSUB  findPnode 

IF  (TestV%  AND  1)=0  THEN  GOSUB  readv  ’  Set  all  potentials  to  Zero,  (file  may  have  junk!) 

'  Do  this  first  because  vpufe  and  other  similar  variables  lost  after  CHAIN  to  Connects 

menu0%=3 :  menul%=l 

IF  (TcstE%  AND  1)=0  THEN  CHAIN  Connects 

’IF  TestE%  =0  OR  TestE%=2  THEN  CHAIN  Connects 

IF  (TestE%  AND  2)=0  THEN  CHAIN  WriteEquS 

TF  TestE%<2  THEN  CHAIN  WriteEquS 

menu0%=0:  menul%=0 

CHAIN  SolveS 

**•••*•••••»  CASE  4  Examine  File  Contents  •••********•• 

CASE  4  '  Examine  Menu  selected 

IF  menu  1% <9  AND  MultiCommand%=0  THEN 
IF  Printf%  THEN  changeCursor  4 

'  IF  MeshFile%=0  THEN  PRINT  "File  Not  Open." :  GOTO  Idle 
IF  MeshFile%=0  THEN 
MultiCommand%=  1 
mcnu0p%=4 
menu  1  p%=menu  1  % 
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menuO%=  1  :  menul%=l :  GOTO  scase 
END  IF 

IF  Nodes9b=0  AND  Elements%=0  THEN  PRINT  "File  Is  Empty." :  GOTO  Idle 
END  IF 

SELECT  CASE  menu  19k 

CASE  1  ’  Examine  Whole  File 
MultiCommand%= 1 
menu0p%=4 
menulp9b=l 

IF  level%=0  THEN  level%=l :  menul9k=4  :  GOTO  scase 
IF  level%=  1  THEN  level%=2 :  menul9k=2 :  GOTO  scase 
IF  level%=2  THEN  level%=3  :  menul9k=3  :  GOTO  scase 
IF  level%=3  THEN  level9b=4  :  menul%=5  :  GOTO  scase 
IF  level%=4  THEN  level%=5  :  menul%=6 :  GOTO  scase 
IF  level%=5  THEN  level%=6 :  menul%=7  :  GOTO  scase 
IF  levcl%=6  THEN  level%=7 :  menul9k=8  :  GOTO  scase 
MultiCommand9k=0 

CASE  2  '  Examine  Node  Positions 
IF  xyF%=0  THEN  GOSUB  readxy 
PRINT  #2,  "Xmin  =";xmin;TAB(20);"Xmax  =";xmax 
PRINT  #2,  "Ymin  =";ymin;TAB(20);"Ymax  =";ymax 
PRINT  #2,  "Node";TAB(10);"X";TAB(26);"Y" 

FOR  i9k=l  TO  Nodes9k 
PRINT  #2,  i%;TAB(9);x(i%);TAB(25);y(i%) 
menu0%=MENU(0) :  IF  menu0%>0  THEN  m select 
NEXT i% 

CASE  3  '  Examine  Node  Potentials 
IF  vF%=0  THEN  GOSUB  readv 
PRINT  #2,  "Vmin  =";vmin;TAB(22);"Vmax  =";vmax 
PRINT  #2.  "Node";TAB(10);"V" 

FOR  i%=l  TO  Nodes% 

PRINT  #2,  i9k;TAB(9);v(i%) 
menu09b=MENU(0) :  IF  menu09k>0  THEN  mselect 
NEXT  i% 

CASE  4  '  Examine  Node  Connections 
IF  nuvcF9k=0  THEN  GOSUB  readnuvc 
PRINT  #2,  "Node";TAB(10);"Nu";TAB(19);"Vc" 

FOR  i%=l  TO  Nodes% 

PRINT  #2,  i9k;TAB(9);Nu9k(i%);TAB(18);Vc(i%) 
menu09k*MENU(0) :  IF  menu0%>0  THEN  mselect 
NEXT  i% 

CASES  ’  Examine  Element  Composuon 
IF  nijklF%=0  THEN  GOSUB  readnijkl 
’  IF  nregF%=0  THEN  GOSUB  Findnreg 
HUNT  #2, 

"Element"  ;TAB(10);"Region";TAB(20);"Ni";TAB(29);"Nj";TAB(38);"Nk";TAB(47);"Nl" 

GET  #  1 ,  regptA 

FOR  i%=l  TO  Elements% 

GET  #1 

Region%=CVI(aa2S) 

PRINT  #2, 

i%;TAB(9)0<egion%;TAB(19).Ni%(i%);TAB(28);Nj%(i%);TAB(37);Nk%(i%);TAB(46);NI%(i9k) 
menu09k=MENU(0) :  IF  menu09b>0  THEN  mselect 
NEXT  i% 

CASE  6  '  Examine  Element  Properties 
IF  mcurdF%=0  THEN  GOSUB  readmcurd 
IF  uF%=0  THEN  GOSUB  readu 
IF  typeF%=0  THEN  GOSUB  FindType 
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1  *  Tr*'« 

"Elemcm";TAB(10);"Matc";TAB(16);"Type”;TAB(22);"Ciird";TAB(37);"Pennr;TAB(52);"Penn2 

FOR  i%=l  TO  Elements% 

PRINT  #2,  i%;TAB(9);mate%(i%);TAB(15);Type%(i%);TAB(21);Curd(i%); 

PRINT  #2,  TAB(36);ul(i%);TAB(51);u2(i%) 
menuO%=  MENU (0) :  IF  menu0%>0  THEN  msclect 
NEXT i% 

CASE  7  ’  Examine  Element  Connections 
IF  nabcdF%=0  THEN  GOSUB  FindNabcd 

PRINT #2,  "Element";TAB(10);"xNa";TAB(19);"xNb";TAB(28);"xNc";TAB(37);"xNd" 

FOR  i%=l  TO  Elements% 

PRINT  #2,  i%;TAB(9);Na%(i%);TAB(18);Nb%(i%);TAB(27);Nc%(i%);TAB(36);Nd%(i%) 
menu0%=MENU(0) :  IF  menu0%>0  THEN  m select 
NEXT i% 

CASE  8  ’  Examine  Solution  Equations 
IF  TestE%<2  THEN  CHAIN  WriteEquS 
GET#l,equpt&  1  StartVequ&  StartBequ& 

StartVequ&=CVL(ca4$) 

StartBequ&=CVL(cb4$) 

'  PRINT  "From  equpufe  =";equpt&;"  Find  StartVequ&  =";StartVequ& 

GET#1,  StartVequ&  ’  inodes%,  blank%,  Slaitlequ& 
inodes%=CVI(ba2$) 

StartDequ&=CVL(bc4$) 

PRINT  #2,  inodes%;"Set  Node  Potentials  starting  at";StartVequ& 

FOR  i%=l  TO  inodes% 

GET  #1 

Vi%=CVI(ba2$) :  cl=CVS(bc4$) 

PRINT #2,  ”V"+MID${STR$(Vi%)>2>+"=";cl 
menu0%=MENU(0) :  IF  menu0%>0  THEN  m select 
NEXT i% 

GET  #1,  StartDequ&  1  inodes%,  blank%,  StartDequ& 
inodes%=CVI(ba2$) 

StartIequ&=CVL(bc4$) 

PRINT  #2,  modes'Jfc;’’ Dependent  Node  Equations  starting  at";StartDequ& 

FOR  i%=l  TO  inodes% 

GET  #1 

Vi%=CVI(ba2$) :  Vn%=CVI(bb2$) :  cl*CVS(bc4$) 

PRINT  #2,  " V”+MIDS(STR$(Vi%)^ )+"=V"+MID$(STRS(Vn%)^); 

IF  cl<0  THEN  PRINT  #2,  cl  ELSE  PRINT  #2,  ’V+MID$(STRS(cl),2) 
menu0%=MENU(0) :  IF  menu0%>0  THEN  m select 
NEXT  i% 

GET  #  1 ,  Startlequ&  ’  inodes%,  blank%,  StartDequ& 
inodes%=CVI(ba2$) 

StartBequ&=CVL(bc4$) 

PRINT  #2,  inodes%;"Independent  Node  Equations  starting  at";StartIequ& 

FOR  i%=l  TO  inodes% 

GET  #1 

Vi%=CVI(ba2$) :  Nline%=CVI(bb2S) :  cl=CVS(bc4S) 

IF  cl=0  THEN  CS=”  ELSE  C$=STRS(cl) 

PRINT  #2,  "V"+MID$(STR$(Vi%)2)+"=("+C$; 

FOR  j%=  1  TO  Nline% 

GET  #1 

Vn%=CVI(ba2S) :  un%=CVI(bb2$) :  kij=CVS(bc4$) 

IF  kij=l  THEN  k$=""  ELSE  k$=STRS(kij) :  IF  kij>=0  THEN  MID$(k$,l,l)="+" 

IF  un%=0  THEN  u$=""  ELSE  u$="u"+MID$(STRS(un%),2) 

IF  k$=""  THEN  uS="+"+uS  ELSE  uS="*"+uS 
v$="V"+MIDS(STRS(Vn%),2) 

IF  u$o"+"  AND  u$o"*"  THEN  v$="*"+v$ 

PRINT  #2,  k$+uS+vS; 
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IF  un%=0  AND  ldj=l  THEN  d$=d$+’l"  ELSE  d$=dS+k$+u$ 

NEXT  j% 

IF  d$=T  THEN  PRINT  #2,  ")"  ELSE  PRINT  #2,  ")/("+d$+T 
menu0%=MENU(0) :  IF  menu0%>0  THEN  m select 
NEXT  i% 

GET#1,  StartBequ&  '  inodes%,  blank%,  Endequ& 
inodes%=CVI(ba2$) 

PRINT  #2,  inodes%;"Flux  Density  Equations  starting  at";StartBequ& 

FOR  i%=l  TO  inodes% 

GET  #1 

nelem%=CVI(aa2$) :  vli%=CVI(ab2$) :  vlj%=CVI(ac2S) :  viv?i_C  vi(ad2$) 

GET  #1 

ci=CVS(ca4$) :  cj=CVS(cb4$) :  ck=-ci-cj 
GET  #1 

bi=CVS(ca4$) :  bj=CVS(cb4$) :  bk=-bi-bj 
PRINT  #2  "Bx- 

E"+MID$(STR$(nelem%),2)+"=";ci;"*V"+MID$(STR$(vli%),2)+"+";cj;"*V”+MID$(STR$(vlj%),2)+"+ 

";ck;"*V"+MID$(STR$(vlk%),2) 

PRINT  #2,  "By- 

E"+MID$(STRS(nelem%),2)+"=";bi;"*V"+MID$(STRS(vli%)^)+"+";bj;"*V',+MIDS(STRS(vlj%)^)+”+ 
"  ;bk; "  *  V"+MIDS(STR$(  v  1  k%),2) 

menu0%=MENU(0) :  IF  menu0%>0  THEN  m select 
NEXT  i% 

CASE  10 
CLOSE  #2 

OPEN  "SCRN:"  FOR  OUTPUT  AS  #2 
Prinlf%=0 

CASE  11 

mne$=FILES$(0,"Disk  File  for  Screen  Dump") 

IF  mne$<>""  THEN 
CLOSE  #2 

OPEN  mneS  FOR  OUTPUT  AS  #2 
Printf%=l 
END  IF 

CASE  12 

'  Cancel  Print,  (goto  Idle) 

MultiCommand%=0 

CASE  ELSE 

PRINT  "Menu  4  Not  Finished  Yet" :  INPUT  x  :  STOP 
END  SELECT 

IF  menul%<9  AND  MultiCommand%=0  THEN 
PRINT  #2,  "Eid  of  File." 

END  IF 


CASE  5  Plot  Material  Outlines,  Mesh,  and  Flux  (constant  potential) 


CASE  5  '  Plot  selected 
IF  MeshFiJe%=0  THEN 
MuItiCommand%=  1 
menu0p%=5 
menu  1  p%=menu  1  % 

menu0%=l  :  menul%=l :  GOTO  scase 
ELSE 

IF  PloiFluxF%=0  THEN  MultiCommand%=0  '  Draw  Material  Outline  with  fluxplot 
END  IF 

IF  Nodes%=0  AND  Elements%=0  THEN  PRINT  "File  Is  Empty." :  GOTO  Idle 
changeCursor  4 
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IF  nabcdF%=0  AND  (menul%=l  OR  menul%=4)  THEN  GOSUB  FindNabcd 
IF  xyF%=0  THEN  GOSUB  readxy 
IF  nijklF%=0  THEN  GOSUB  readnijkl 


SELECT  CASE  menul% 

CASE  1  1  Draw  Material  Outline 

IF  mcurdF%=0  THEN  GOSUB  readmcurd 
IF  typeF%=0  THEN  GOSUB  FindType 
frame  xmin,xmax,ymin,ymax,-.9 
CLS 

PRINT  xmin,xmax,ymin,ymax 
mate%(0)=-l 

FOR  i%=l  TO  Elements% 
ml%=mate%(i%) 

IF  mate%(Na%(i%))om  1  %  THEN  draw  x(Ni%(i%)),  y(Ni%(i%)),  x(Nj%(i%)),  y(Nj%(i%)) 
IF  mate%(Nb%(i%))oml%  THEN  draw  x(Nj%(i%)).  y(Nj%(i%)),  x(Nk%(i%)),  y(Nk%(i%)) 
IF  mate%(Nc%(i%))oml%  THEN  draw  x(Nk%(i%)),  y(Nk%(i%)),  x(Nl%(i%)),  y(Nl%(i%)) 
IF  mate%(Nd%(i%))oml%  AND  Type%(i%)>=0  THEN  draw  x(Nl%(i%)),  y(Nl%(i%)), 
x(Ni%(i%)),  y(Ni%(i%)) 

menu0%=MENU(0) :  IF  menu0%>0  THEN  m select 
NEXT  i% 

CASE  2  ’  Draw  Mesh 

IF  typeF%=0  THEN  GOSUB  FindType 
frame  xmin,xmax,ymin,ymax,--9 
CLS 

PRINT  xmin,xmax,ymin,ymax 
FOR  i%=l  TO  Elements% 

nl%=Ni%(i%):  n2%=Nj%(i%) :  n3%=Nk%(i%):  n4%=Nl%(i%) 

SELECT  CASE  Type%(i%) 

CASEO 

draw  x(nl%),  y(nl%),  x(n2%),  y(n2%) 
drawto  x(n3%),  y(n3%) 
drawto  x(n4%),  y(n4%) 
drawto  x(nl%),  y(nl%) 

CASE  1 

draw  x(n4%),  y(n4%),  x(nl%),  y(nl%) 
drawto  x(n2%),  y(n2%) 
drawto  x(n3%),  y(n3%) 
drawto  x(n4%),  y(n4%) 
drawto  x(n2%),  y(n2%) 

CASE  2 

draw  x(nl%),  y(nl%),  x(n2%),  y(n2%) 
drawto  x(n3%),  y(n3%) 
drawto  x(n4%),  y(n4%) 
drawto  x(n  1  %),  y(n  1  %) 
drawto  x(n3%),  y(n3%) 

END  SELECT 

menu0%=MENU(0) :  IF  menu0%>0  THEN  mselect 
NEXT i% 

CASE  3  '  Draw  Region 

IF  nregF%=0  THEN  GOSUB  Findnreg 
PRINT  "Draw  Region  not  yet  finished" 

CASE  4  ’  Draw  Flux  Plot 

IF  vF%=0  THEN  GOSUB  readv 

IF  vmax-vmin=0  THEN  PRINT  "No  Potentials  Defined" :  GOTO  Idle 
IF  PlotFluxF%=0  THEN 
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INPUT  "Flux  lines  to  plot",Nline% 

IF  Nline%>200  THEN  Nline%=200 

IF  Nline%=0  THEN  Nline%=10 

PlotFluxF%=l 

MultiCommand%= 1 

menu0p%=5 

menulp%=4 

menu0%=5  :  menul%=l :  GOTO  scase 
ra 

PlotFluxF%=0 
MultiCommand%=0 
END  IF 

FOR  i%=l  TO  Nline% 
plot(i%>=vmin+{vmax-vmin)*(i%-.5)/Nline% 

NEXT  i% 

FOR  i%=l  .  0  Elements% 

nl%=Ni%(i%) :  n2%*Nj%(i%) :  n3%=Nk%(i%) :  n4%=Nl%(i%) 

vl=v(nl%) :  v2=v(n2%) :  v3=v(n3%) :  v4=v(n4%) 

n5%=nl%  :  n6%=n2%  :  n7%=n3%  :  n8%=n4% 

v5=vl :  v6=v2  :  v7=v3  :  v8=v4 

IF  vl>v6 THEN  SWAP  vl,v6:  SWAPnl%,n6% 

IF  v2>v7  THEN  SWAP  v2,  v7  :  SWAP  n2%,  n7% 

IF  v3>v8  THEN  SWAP  v3,v8:  SWAPn3%,n8% 

IF  v4>v5  THEN  SWAP  v4.  v5  :  SWAP  n4%,  n5% 

FOR  j%=l  TONline% 
test%=0 
v=plot(j%) 

IF  v>v  1  AND  v<v6  THEN 
test%=test%+l 
IF  test%=l  THEN 
Inter  vl,x(nl%),v6,x(n6%),v.xl 
Inter  vl,y(nl%),v6,y(n6%),v,yl 
ELSE 

Inter  v  1  ,x(n  1  %),v6pc(n6%)  ,v,x2 
Inter  vl,y(nl%),v6.y(n6%),v,y2 
draw  xl,ylpi2,y2 
test%=0 
END  IF 
END  IF 

IF  v>v2  AND  v<v7  THEN 
test%=test%+l 
EFtest%=l  THEN 
Inter  v2vx(n2%)  ,v7  ,x(n7  % )  ,v  ,x  1 
Inter  v2,y(n2%),v7,y(n7%),v,yl 
ELSE 

Inter  v2,x(n2%),v7,x(n7%),v,x2 
Inter  v2,y(n2%),v7,y(n7%),v,y2 
drawxl,yl,x2,y2 
iest%=0 
END  IF 
END  IF 

IF  v>v3  AND  v<v8  THEN 
test%=test%+l 
IF  test%=l  THEN 
Inter  v3,x(n3%),v8,x(n8%),v,xl 
Inter  v3,y(n3%),v8,y(n8%),v.yl 
ELSE 

Inter  v3,x(n3%),v8,x(n8%),v,x2 
Inter  v3,y(n3%),v8,y(n8%),v,y2 
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draw  xl,yl,x2,y2 
test%=0 
END  IF 
END  IF 


IF  v>v4  AND  vcv5  THEN 
test%=test%+l 
IF  lest%=l  THEN 
Inter  v4,x(n4%),v5,x(n5%),v.xl 
Inter  v4,y(n4%),v5,y(n5%),v,yl 
ELSE 

Inter  v4,x(n4%),v5rx(n5%),v,x2 
Inter  v4,y(n4%),v5,y(n5%),v,y2 
draw  xl,yl,x2,y2 
test%=0 
END  IF 
END  IF 

NEXTj% 

NEXT  i% 

END  SELECT 


***••»»*****  CASE  6  Integrate  H  •  dl  **•*•*••***•* 

CASE  6  '  Integrate  selected 

IF  MeshFilc%=0  THEN 
MultiCommand  %=  1 
menu0p%=6 
menulp%=l 

menuO%=l :  menul%=l :  GOTOscase 
END  IF 
changeCursor  4 

IF  vF%=0  THEN  GOSUB  rcadv 

IF  vmax-vmin=0  THEN  PRINT  "No  Potentials  Defined"  :  GOTO  Idle 

'  Trashing  info  in  Nu%0,  VcO,  plotO 
IF  xyF%=0  THEN  GOSUB  readxy 
IF  nijklF%=0  THEN  GOSUB  readnijkl 
IF  typeF%=0  THEN  GOSUB  FindType 
IF  uF%=0  THEN  GOSUB  readu 

Hagain:  WINDOW  OUTPUT  1 
IF  PlotFluxF%=0  THEN 
PlotFluxF%=l 
MultiCommand%= 1 
menu0p%=6 
menulp%=l 

menu0%=5  :  menul%=l :  GOTO  scase 
ELSE 

PlotFluxF%=0 
END  IF 

MultiCommand%=0 

INITCURSOR 

Hdl#xO# 

WINDOW  2„(35030H500^00),4  '  512340  bottom  right 
CLS 

PRINT  "Integrate  H*dl" 

PRINT  "(xl,yl)to(x2,y2)" 

INPUT  "xl";x2 
INPUT  "yl";y2 
m point:  xl*x2  :  yl=y2 
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INPUT  "x2";x2 
INPUT  ”y2";y2 

WINDOW  OUTPUT  1  '  Don't  cover  window  2 

CALL  defnodecond  '  Define  position  definition  in  Nu%  for  each  node 
IF  menu0%>0  THEN  m  select 
FOR  n%=l  TO  2 
begin  %=1 
WHILE  begin  %>0 

getdlelem  n%,begin%,elem%,xhl1ynl,xh2,yh2, factor 
EF  elem%>0  THEN 

draw  xhl,yhl,xh2,yh2 
xli=x(nll%):  yli=y(nll%) 
xlj=x(nl2%) :  ylj=y(nl2%) 
xlk=x(nl3%) :  ylk=y(nl3%) 

al4=l/((xli-xlj)*(yli-ylk)+(yii-ylj)*(xlk-xli))  '  l/(2*Area)  positive  direction. 

bli=(ylj-ylk)*al4 :  bli=(ylk-yli)*al4 :  blk=-bli-blj 

cli=(xlk-xlj)*al4 :  clj=(xli-xlk)*al4 

Bx=(v(n  1 1  %)-  v(n  1 3  %))*c  1  i+(  v(n  1 2%)-v(n  1 3  %))*c  1  j 

By=(v(nl3%)-v(n  1 1  %))*bl  i+(v(n  13%)-v(n  12%))*bl  j 

IF  n%=l  THEN 

Hx=Bx*ul(ekm%)  ’  ul  is  \t\x  (inverse  of  relative  pennubility 
Hy=By*ul(elem%) 

ELSE 

Hx=Bx*u2(elem%) 

Hy=By*  u2(eletn%) 

END  IF 

Hdl#=Hdl#+(Hx*(xh2-xhl)+Hy*(yh2-yhl))*factor 
•PRINT  "Bx  =";Bx 
-PRINT  "By  =";By 
END  IF 

begin%=elem%+l 
WEND 
NEXT  n% 

WINDOW  OUTPUT  2 
PRINT  "H*dl=";CSNG(Hdl#AiO) 

INPUT  "Continue  (y  ji  or  r)”;x$ 

IF  UCASES(MID$(x$,  1 , 1  ))="R  ”  THEN  Hagain  ’  Start  over  with  another  integration 
IF  UC ASE$(MID$(x$,  1 , 1  ))o"N"  THEN  mpoint  '  Continue  integration  to  another  point 
MultiCotnmand%=0 
WINDOW  1 


CASE  ELSE 


CASE  ELSE 

PRINT  "Opps !!!  Something  wierd  has  happened.  Hit  RETURN  and  proceed  with  caution." 
MulnCommand%=0 
INPUT  x 
END  SELECT 

Idle:  IF  MultiCommand%  THEN 
menuO%=menuOp% 
menu  1  %smenu  lp% 

ELSE 

menu0%=0 

INTTCURSOR 

MENU  1,1,1, "Open"  '  Deselect  Menu  Items  —  Idle  loop 
END  IF 

WEND 
'  Idle  Stuff 

IF  MeshFile%  THEN  '  Read  File  while  waiting  for  user  input. 

IF  xyF%=0  THEN  GOSUB  readxy 
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IF  nijklF%=0  THEN  GOSUB  readnijkl 
IF  vF%=0  THEN  GOSUB  readv 
IF  mcurdF%=0  THEN  GOSUB  readmcurd 
IF  uF%=0  THEN  GOSUB  readu 
IF  typcF%=0  THEN  GOSUB  FindType 

IF  nabcdF%=0  AND  (TestE%=l  OR  TestE%=3)  THEN  GOSUB  FindNabcd 
IF  nuvcF%=0  THEN  GOSUB  readnuvc 
IF  nregF%=0  THEN  GOSUB  Findnreg 
END  IF 
GOTO  loop 

'  ###################  subs  mntimMttttmttHttttttM 

readnuvc: '  Node  Nu  Vc" 

•  PRINT  "Reading  Nu,  Vc" 

WHILE  inu  vc%<Nodes  1  % 

GET  #1,  inuvc%+l 
Nu%(inuvc%>=CVI(ba2$) 

Vc(inuvc%)=CVS(bc4S) 
inu  vc%=inuvc%+ 1 

menu0%=MENU(0) :  IF  menu0%>0  THEN  m select 
WEND 
inuvc%=l 
nuvcF%=l 
RETURN 

rcadxy: '  Node  X  Y" 

PRINT  "Reading  X,  Y" 

'  xypt&=l+Nodes% 

IF  ixy%=l  THEN 
GET  #1,  l+xypt& 

xmin=CVS(ca4$) :  xmax=CVS(ca4S) 
ymin=CVS(cb4$) :  ymax=CVS(cb4$) 

END  IF 

WHILE  ixy%<Nodesl% 

GET  #1,  ixy%+xypt& 

x=CVS(ca4S) 

y=CVS(cb4$) 

IF  x<xmin  THEN  xmin=x 

IF  y<ymin  THEN  ymin=y 

IF  x>xmax  THEN  xmax=x 

IF  y>ymax  THEN  ymax=y 

x(ixy%>=x 

y(ixy%)=y 

ixy%=ixy%+l 

menu0%=MENU(0) :  IF  menu0%>0  THEN  m select 
WEND 

LSET  ca4S=MKS$(xmin) :  LSET  cb4S=MKS$(xmax) 

PUT  #l,equpt&-fl 

LSET  ca4S=MKS$(ymin) :  LSET  cb4S=MKSS(ymax) 

PUT  #1,  equpt&+2 
ixy%=l 
xyF%=l 
RETURN 

readv: 

1  PRINT  "Reading  V" 

’  vpt&=l+2*Nodes% 

IF  (TestV%  AND  1)  THEN 

IF  iv%=l  THEN  GET  #l,l+vpt&  :  vmin=CVS(ca4S) :  vmax=CVS(ca4$) 
WHILE  iv%<Nodesl% 

GET  #1,  iv%+vpt& 
v=CVS(ca4$) 

IF  v<vmin  THEN  vmin=v 
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IF  v>vmax  THEN  vmax=v 
v(iv%)=v :  iv%=iv%+l 

menuO%=MENU(0) :  IF  menu0%>0  THEN  mselect 
WEND 
ELSE 

WHILE  iv%<Nodesl% 
v(iv%)=0  '  Stuff  with  0  Potentials 
LSET  ca4$=MKS$(0) 

PUT  #1,  iv%+vpt& 
iv%=iv^fc+l 

menu0%=MENU(0) :  IF  menu0%>0  THEN  mselect 
WEND 

'  HUNT  "Potentials  initialized  to  zero" 
vmin=0:  vmax=0 
TestV%=TestV  %  OR  1 
GET  #1,1 

LSET  ac2£=MKIS(TestV%) 

PUT  #1, 1 
END  IF 

LSET  ca4S=MKSS(vmin) :  LSET  cb4S=MKSS(vmax) 

PUT  #1,  equpt&+3 

iv%=l 

vF%=l 

RETURN 

readnijkl: '  Element  Ni  Nj  Nk  Nl" 

'  HUNT  "Reading  Ni,  Nj,  Nk,  Nl" 

’  nijklpt&=  l+3*Nodes% 

WHILE  inijkl%<Elementsl% 

GET  #1,  inijkl%+nijklpt& 

Ni%(inijkl%>=CVI(aa2$) 

Nj%(inijkl%)=CVI(ab2$) 

Nk%(inijkl%)=CVI(ac2S) 

Nl%(inijkl%>=CVI(ad2$) 

inijkl%=inijU%+l 

menu0%=MENU(0) :  IF  menu0%>0  THEN  mselect 
WEND 
inijkl%sl 
nijklF%=l 
RETURN 


readme urd: '  Element  Mate  Curd 
'  PRINT  "Reading  Mate,  Curd" 

'  mcuidpt&=l+3*Nodes%+Elements% 

WHILE  imcurd%<Elementsl% 

GET  #1,  imcurd%+mcurdpt& 
mate%(imcurd%)=CVI(ba2S) 
Curd(imcurd%)=CVS(bc4$) 
imcurd%=imcurd%+ 1 

menu0%=MENU(0) :  IF  menu0%>0  THEN  mselect 
WEND 
imcurd%*l 
mcurdF%=l 
RETURN 

readu: '  Element  Perm 
HUNT  "Reading  u I,  u2" 

'  upt&=l+3*Nodes%+2*Elements% 

WHILE  iu%<Elementsl% 

GET  #1,  iu%+upt& 
u  1  (iu%)=CV  S(ca4$) 
u2(iu%>*CV  S(cb4$) 
iu%*iu%+l 

menuO%=MENU(0) :  IF  menu0%>0  THEN  mselect 
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WEND 
iu%=  1 
uF%=l 
RETURN 

Findnreg:  *  Element  Region 
'  PRINT  "Finding  number  of  regions" 

'  regpt&=l+3*nodes9b+3*elements% 

WHILE  ireg%<Elementsl% 

GET  #1,  neg%+regpt& 

Region%=CVI(aa2$) 

IF  Region%>nreg%  THEN  nreg$>=Region% 
ireg%=ircg%+ 1 

menuO%=MENU(0) :  IF  menu0%>0  THEN  m select 
WEND 
ircg%=l 
nregF%=l 
RETURN 

FindNabcd: ' 

'  PRINT  "Creating  Triagonal  element  types  from  square  types  where  appropriate" 
'  nabcdpt&=l+3*Nodes%+4*Elements% 

WHILE  inabcd%<Elementsl% 

GET  #1,  inabcd%+nabcdpt& 

Na%(inabcd%)=CVI(aa2$) 

Nb%(inabcd%)=CVI(ab2$) 

Nc%(inabcd%)=CVI(ac2$) 

Nd%(inabcd%)=CVI(ad2$) 
inabcd%=inabcd%+ 1 

menu0%=MENU(0)  :  IF  menu0%>0  THEN  m select 
WEND 
inabcd%=l 
nabcdF%=l 
RETURN 

FindType: 

’  PRINT  "Splitting  Elements  into  triangular  where  appropriate" 

IF  xyF%=0  THEN  GOSUB  readxy 
IF  nijklF%=0  THEN  GOSUB  readnijkl 
WHILE  itype%<Elementsl% 

GET  #1,  itype%+typept& 

Type%=CVI(bb2$) 

IT  TytyflLg/)  TUpKT 

nl%=Ni%(i%):  n2%=Nj%(i%)  :  n3%=Nk%(i%) :  n4%=Nl%(i%) 

xi=x(nl%) :  yi=y(nl%) 

xj=x(n2%):  yj=y(n2%) 

xk=x(n3%) :  yk=y(n3%) 

xl=x(n4%) :  yl=y(n4%) 

d  1  »(xk-xi)A2+<yk-yi)A2 

d2»(xl-xj)A2+(yl-yj)A2 

Type%*l 

IF  dl<d2*.9999  THEN  Type%=2 
LSET  bb2S*MKIS(Type%) 

PUT  #1,  itype%+typept& 

'  PRINT  i%,Type% 

END  IF 

Type%(itype%)=Type% 
itype%=itype<&+ 1 

menu0%=MENU(0) :  IF  menu0%>0  THEN  m select 
WEND 
itype%=l 
typeF%=l 
RETURN 
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flndPnode: 

'SHARED  xyF%,Pnode%,StartVequ&,equpt&,TestV% 

IF  xyF%=0  THEN  GOSUB  readxy 
PRINT  .07„08 

INPUT  "Coordinates  of  Node  to  monitor  during  solve:  x,  y  \xn,yn 
dmin=(x(l)-xn)A2+(y(l)-yn)A2 :  Pnode%=l 
FOR  i%=2  TO  Nodes% 
d=(x(i%)-xn)A2+{y(i%)-yn)A2 
IF  d<dmin  THEN  dmin=d :  Pnode%=i% 

NEXT  i% 

PRINT  "Will  monitor  Node";Pnode% 

StartVequ&=equpt&+4  '  ***  Potential  equations  start  here 

GET  #1,  StartVequ& 

LSET  bb2$=MKIS(Pnode%)  '  Node  to  prim  out 
PUT  #1,  StartVequ& 

TestV%*TestV%  OR  2 
GET  #1,1 

LSET  ac2$=MKI$(Test  V  %) 

PUT  #1,1 
RETURN 


|  GRAPHIC  ROUTINES  FOR  DEFAULT  WINDOW 

'  Define  graphic  coordinates  of  default  window  (quadrant  1  coordinates). 

’  xl,x2,yl,y2  =  coordinates,  xl<x2,  yl<y2 

ABS(windowfraction)  =  is  %  of  window  that  the  coordinates  xl,x2,yl,y2 
represent  (0.9  is  a  good  number  to  use) 

'  NOTE:  windowfractioncl  forces  linear  x-y  scaling 

SUB  frame(xl,x2,yl,y2.windowfraction)  STATIC 
SHARED  xmuldplierscale.xshiftscale.ymultiplierscale.yshiftscalc 
test=  1  -SGN(  windowfraciion) 
windowfracuon=ABS(windowfraction) 
wx%=WINDOW(2)-l  '  wirdow  width 
Hy%=WINDOW(3)-l  '  window  height 
xl0sINT((l-windowfraction)*wx%/2)  ’  SYSTEM(5)=screen  width 
y  10=INT((  1  -wmdowfraction)*  Hy %/2)  '  SYSTEM(6)=screen  height 
x20swx%-xl0 
y20=Hy%-y!0 

xmultiplierscale=(x20-x  1 0V(x2-x  1 ) 
ymuluplierscale=(y  10-y20)/(y2-y  1 ) 
xmult=ABS(xmultiplierscale) 
ymult=ABS(ymultiplierscale) 

IF  test  AND  xmuloymult  THEN  xmultipliersca)e=ymult*SGN(xmultipliersca)e) 
IF  test  AND  ymult>xmult  THEN  ymultiplierscale=xmult*SGN(ymul  tiplierscale) 
xshiftscalesx  1 0-xmul  tiplierscale*  x  1 
yshiftscale*y20-ymul  tiplierscale*  y  1 
END  SUB 


'  Move  pen  to  position  speeded  by  the  frame  subroutine  coordinate  system. 
SUB  pento(xl,yl)  STATIC 

SHARED  xmultiplierscale^sh  if tscale.ymul  tiplierscale, yshiftscale 
MOVETO  x  1  *xmultiplierscale+xshiftscale,y  1  *ymultiplierscale+yshiftscale 
END  SUB 


'  Draw  pen  to  position  speeded  by  the  frame  subroutine  coordinate  system. 

i 

SUB  drawto(x2,y2)  STATIC 
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SHARED  xmultiplierscale.xshiftscale,ymultiplierscale,yshiftscale 
LINETO  x2*xmultiplierscale+xshiftscale,y2*  ymultiplierscale+yshiftscale 
END  SUB 


'  Draw  a  line  as  spec  if ed  by  the  frame  subroutine  coordinate  system. 

t 

SUB  draw(xl,yl,x2,y2)  STATIC 
SHARED  xmuldpUerscale^ishiftscale.ymultiplierscak.yshiftscale 
MO  VETO  x  1  *xmultiplierscale+xshiftscale,y  1  *  ymultiplierscale+yshiftscale 
LINETO  x2*xmultiplierscale+xshiftscale,y2*ymultiplieiscale+yshiftscale 
END  SUB 

'  Compare  two  real  values  to  roughly  the  first  6  digits 
'  1%=0  aob 
•  1%=-1  a=b 

SUB  comparc(a,bJ%)  STATIC 

a&=VARPTR(a) :  POKEL  a&,  PEEKL(aA)  +  128  :  POKE  a&+3, 0 
a&=VARPTR(b) :  POKEL  a&,  PEEKL(a&)  +  128 :  POKE  a&+3, 0 
l%=(a=b) 

END  SUB 

'  Inter  (xl,yl,x2,y2,x,y) 

'  Linear  interpolation  2 

'  Find  y  given  x  and  two  data  points  (xl.yl),  (x2,y2) 

SUB  Inter(x  1  ,y  1  ,x2,y2,x,y)  STATIC 
y=(x-x  1  )/(x2-x  1  )•  (y2-y  1  >+y  1 
END  SUB 

'  Intersect  (xal ,yal ,xa2,ya2,xbl ,ybl ,xb2,yb2,x,y,test%) 

'  Given  two  lines,  one  passing  through  (xal  ,yal),  (xa2,ya2) 

’  and  the  other  passing  through  (xbl.ybl),  (xb2,yb2), 

'  Find  point  of  intersection,  fx.y) 

'  Lines  must  have  non-zero  length,  but  may  be  horizontal  or  vertical 
'  Test%=0  -  No  intersection 
'  Test%=l  -  One  point  of  intersection 
'  Test%=2  -  Infinite  points  of  intersection 
'  Modification:  (x.y)  must  be  on  line  segment  'a'  for  test%>0 

SUB  Intersect(xal,yal,xa2,ya2.xbl,ybl,xb2,yb2,x,y,test$))  STATIC 
dxa=xa2-xal 
dxb=xb2-xbl 
dya=ya2-yal 
dyb=yb2-ybl 
denom=dxa*dyb-dxb*dya 
IF  denom*OTHEN 
testl=dyb*(xal-xbl) 
test2=dxb*(ya  I -y  b  I ) 

a*=VARPTR(testl) :  POKEL  a&,  PEEKL(a&)  +  128  :  POKE  a&+3, 0 
a&=VARPTR(test2) :  POKEL  a&,  PEEKL(a&)  +  128  :  POKE  a&+3, 0 
IF  testl=tesl2  THEN  test9b=2 :  EXIT  SUB 
test%=0 :  EXIT  SUB 
END  IF 
test%=l 

y=(dya*dyb*(xbl-xal>+yal*dxa*dyb-ybl*dya*dxb)/denom  '  point  of  intersection 

x*(dxa*dxb*(ya  1  -yb  1  )-dya*dxb*xal  +dyb*dxa*xb  1  Vdenom 

xl*x*.9999 

x2=x*  1.0001  '••• 

yl=y*.9999  ’*** 

y2*y*  1.0001  '•••toend 

IF  (x2<xal  AND  x2<xa2)  OR  (xl>xal  AND  xl>xa2)  OR  (y2<yal  AND  y2<ya2)  OR  (yl>yal  AND 
yl>ya2)  THEN  test%=0 :  EXIT  SUB 
END  SUB 
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SUB  defnodecond  STATIC 
SHARED  xl.ylA2.y2 

SHAREDnuvcF%^nuvc%,dx,dytNodes%,xyF%fx0.yOJ^u%0^nenuO% 

'  Trashing  info  in  Nu%0 
nuvcF%=0 :  inuvc%=l 
'  Assumes  xyF%=l 

IF  xyF%=0  THEN  PRINT  "ERROR  in  sub  defnodecond:" :  PRINT  "xO  and  yO  not  defined" 
STOP 
dx»x2-xl 
dy=y2-yl 
Nu%(0)=0 

FOR  i%=l  TO  Nodes% 
x=x(i%) 

y=y(i%) 

testl=dy*(x-xl)+(yl-y*1.0001)*dx 
test2=dy*(x-x  l)+(y  1  -y*  .99999)*dx 
IF  testl<0  AND  test2<0  THEN 
'  Point  A  is  to  the  tight  of  a  line  going  from  Point  1  to  Point  2 
test%=l 

ELSEIF  testl>0  AND  tcst2>0  THEN 
'  Point  A  is  to  the  left  of  a  line  going  from  Point  1  to  Point  2 
test%=4 
ELSE 

'  Point  A  is  on  the  line  going  from  Point  1  to  Point  2 
test%=2 
END  IF 

xlt=xl :  x2t=x2 :  ylt=y  1 :  y2t=y2 
'  Round  off  so  real  comparison  works 
a&=VARPTR(x)  '  get  variable  pointer 

POKEL  a&,  PEEKL(a&)  +  128  '  equivent  to  INT(a+.5) 

POKE  a&+3, 0  '  truncate  last  8  bits 

a&=VARPTR(xlt) :  POKEL  a&,  PEEKL(a&)  +  128  :  POKE  a&+3. 0 
a&= V ARPTR(x2t) :  POKEL  a&.  PEEKL(a&)  +  128 :  POKE  a&+3, 0 
a&=VARPTR(y)  :  POKEL  a&,  PEEKL(a&)  +  128 :  POKE  a&+3, 0 
a&=VARPTR(ylt) :  POKEL  a&,  PEEKL(a&)  +  128  :  POKE  a&+3, 0 
a&= V ARPTR(y 2t) :  POKEL  a&,  PEEKL(a&)  +  128 :  POKE  a&+3, 0 
IF  x<=xlt  AND  x<=x2t  THEN 
IF  x=xlt  AND  x=x2t  THEN 
test  %=  test  %  OR  16 
ELSE 

test%=test%  OR  8 
END  IF 

ELSEIF  x>=xlt  AND  x>=x2t  THEN 
test%=test%  OR  32 
ELSE 

test%=test%  OR  16 
END  IF 

IF  y<=ylt  AND  y<=y2t  THEN 
IF  y=ylt  AND  y=y2l  THEN 
test%=test%  OR  128 
ELSE 

test%=test%  OR  64 
END  IF 

ELSEIF  y>=y  1 1  AND  y>=y2t  THEN 
test%=test%  OR  256 
ELSE 

test%=test%  OR  128 
END  IF 

Nu%(i%)=test% 

*  PRINT  i%;test<fc;test%  AND  448;test%  AND  56;test%  AND  7 
menu0%=MENU(0) :  IF  menu0%>0  THEN  WINDOW  1 :  EXIT  SUB 
NEXT  i% 
dx=ABS(dx) 


:  INPUT  x 
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dy=ABS(dy) 
END  SUB 


SUB  getdlelem(n%  .begin  %  ,elem%  ,xh  1  ,yh  1  ,xh2  ,yh2  factor)  STATIC 
'  Assumes  this  routine  will  be  used  directly  after  sub  'defnodecond' 

'  Variables  xl,yl^2ty2,dx,dy^u%0  come  from  defnodecond 
'  Info  in  Nu%0  and  VcQ  is  trashed. 

'  n%=l  or  2.  selects  which  triangle  of  pair  to  investigate 
'  start  seaching  elements  at  begin%  for  intersection  with  line  xl,yl-x2,y2 
'  elem%  is  first  element  found  with  intersection 
'  if  elem%=0  then  search  is  finished 
'  xhl,yhl,xh2,yh2  are  coordinates  of  the  intersection 
'  normally  factors  1. 

'  factors  5  means  the  intersection  is  shared  by  two  elements 
’  nll%,nl2%,nl3%,type%  '  (saves  having  to  look  it  up  later) 

SHARED  Elemenu%fNi%0^j%0^<*>0Jfl%0.Type%0.Vc0.plot0.x0.y0 
SHARED  xl.yl  ,x2,y2,dx,dyfNu%0.memi()% 

SHARED  nl  l%,nl2%,nl3%,Type%  '  (saves  having  to  look  it  up  later) 

FOR  i%=begin%  TO  Elements  % 

nl%=Ni%(i%):  n2%=Nj%(i%):  n3%=Nk%(i%) :  n4%=Nl%(i%) 

Type%=Type%(i%) 

SELECT  CASE  Type% 

CASE  1 

IF n%=l  THEN  nll%=nl%  :  nl2%=n2% :  nl3%=n4%  ELSE  nil %=n2%  :  nl2%=n3%: 
nl3%=n4% 

CASE  2 

IF  n%=l  THEN  nll%=nl%  :  nl2%=n2% :  nl3%=n3%  ELSE  nll%=nl%  :  nl2%=n3%  : 
nl3%=n4% 

CASE  ELSE  'triangle 

IF n%=l  THEN  nll%=nl% :  nl2%=n2% :  nl3%=n3%ELSEnll%=0:  nl2%=0:  nl3%=0 
END  SELECT 

nunll%=Nu%(nll%) :  nunl2%=Nu%(n!2%) :  nunl3%=Nu%(nl3%) 
testl%=nunl  1%  AND  nunl2%  AND  nunl3% 

’  Eliminate  elements  outside  box  volume 

IF  testl%  AND  360 THEN  nextelem  '  8  and  32  and 64  and  256 

test2%=((nunll%  AND  2)=0)  +  ((nunl2%  AND  2)=0)  +  ((nunl3%  AND  2)=0) 

IF  test2%=-3  THEN  '  No  nodes  on  line 
'  Eliminate  elements  above  or  below  line 
IF  testl%  AND  5  THEN  nextelem  '  1  and  4 
factors  1 

ELSEIF  test2%=-2  THEN  '  One  node  on  line 
'  Eliminate  elements  above  or  below  line 
IF  (nunl  1%  AND  2)  THEN 
BF  nun  12%  AND  nunl3%  AND  5  THEN  nextelem 
END  IF 

IF  (nun  12%  AND  2)  THEN 
IF  nunl  1%  AND  nun  13%  AND  5  THEN  nextelem 
END  IF 

IF  (nun  13%  AND  2)  THEN 
IF  nunl  1%  AND  nunl2%  AND  5  THEN  nextelem 
END  IF 
factors  1 

ELSE  '  Two  nodes  on  line 
IF  test2%=0  THEN  nextelem  '  Bad  element 

factors. 5 

END  IF 

'  Element  has  survived  basic  tests--  now  lode  further 
xal=x(nll%):  yal=y(nll%) 
xa2sx(nl2%) :  ya2=y(nl2%) 
xa3=x(nl3%):  ya3=y(nl3%) 

Intersect  xa  1  ,ya  1  ,xa2,ya2,x  1  ,y  1  ,x2,y2,xh  1  ,yh  1  .test  1  % 

Intersect  xa2,ya2,xa3,ya3,xl,yl.x2,y2,xh2,yh2,test2% 

Intersect  xa3  ,ya3  ^ta  1  ,yal  ,x  1  ,y  1  tx2,y2,xh3,yh3,test3% 
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test%sO 

IF  tesil%=l  THEN 
test%=test%+l 

Vc(tesx%)=xh'  :  Vc(test%+8)=yhl 
ELSEIF  tcstl%=2  THEN 
test%=test%+l 

Vc(test%)=xal :  Vc(test%+8)=yal 
test%=tcst%+ 1 

Vc(test%)=xa2 :  Vc(test%+8)=ya2 
ELSE 
END  IF 

IF  test2%=l  THEN 
test%*test%+l 

Vc(test%)=xh2 :  Vc(test%+8)=yh2 
ELSEIF  test2%=2  THEN 
test%=test%+l 

Vc(test%>=xa2 :  Vc(test%+8)=ya2 
test%=test%+l 

Vc(test%)=xa3 :  Vc(test%+8)=ya3 
ELSE 
END  IF 

IF  tcsi3%=l  THEN 
test%=test%+l 

Vc(test%)=xh3 :  Vc(test%+8)=yh3 
ELSEIF  test3%=2  THEN 
test%=test%+l 

Vc(test%)=xa3 :  Vc(test%+8)=ya3 
icst%=tesl%+l 

Vc(test%)=xal :  Vc(test%+8)=yal 
ELSE 
END  IF 

IF  iest%=l  THEN  nextclem 
IFdxxiy  THEN 

’  Make  comparisons  based  on  x  values 
xhmin=Vc(l)  :  nmin%=l 
xhmax=xhmin :  nmax%=l 
FOR  j%=2  TO  test% 
x=Vc(j%) 

IF  x<xhmin  THEN  xhmin=x  :  nmin%=j% 

IF  x>xhmax  THEN  xhmax=x  :  nmax%=j% 

NEXT  j% 

compare  xhmin.xhmax.1% 

IF  1%  THEN  nextclem  '  length  of  line  segment  is  zero 

IF  (xhmax<=xl  AND  xhmax<=x2)  OR  (xhmin>=xl  AND  xhmin>=x2)  THEN  nextclem 
plot(l)=Vc(nmin%+8) :  plot(2)=Vc(nmax%+8) 

Vc(l)=xhmin :  Vc(2)=xhmax 
Vc(3)=xl :  pIot(3)=yl 
Vc(4)=x2 :  plot(4>sy2 
FOR  j%=1  TO  3 
FOR  k%=j%+]  TO  4 
IF  Vc(j%)>Vc(k%)  THEN 
x*Vc(k%):  p=plot(k%) 

Vc(k%)=Vc(j%) :  plot(k%)=plot(j%) 

ND^^x  ' 

NEXT k% 

NEXT  j% 

xhl=Vc(2) :  yhl=plot(2) 
xh2=Vc(3) :  yh2*plot(3) 

IF  SGN(x2-x  1  )oSGN(xh2-xh  1 )  THEN 
'  Make  sure  segment  is  going  in  same  direction  as  line  of  integration! 

SWAP  xhl,xh2 
SWAP  yhl,  yh2 
END  IF 
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ELSE 

‘  Make  comparisons  based  on  y  values 
xhmin=Vc(9)  :  nmin%=l 
xhmax=xhmin :  nmax%=l 
FOR  j%=2  TO  test% 
x=Vc(j%+8) 

IF  x<xhmin  THEN  xhmin=x  :  nmin%=j% 

IF  x>xhmax  THEN  xhmax=x :  nmax%=j% 

NEXTj% 

compare  xhmin,xhinax,l% 

IF  1%  THEN  nextelem  ’  length  of  line  segment  is  zero 

IF  (xhmax<=yl  AND  xhmax<=y2)  OR  (xhmin>=yl  AND  xhmin>=y2)  THEN  nextelem 
plot(l)=Vc(nmin%) :  plot{2)=Vc(nmax%) 

Vc(l)=xhmin :  Vc(2)=xhmax 
Vc(3)=yl :  plot(3)=xl 
Vc(4)=y2 :  plot(4)=x2 

■PRINT  ”%";Vc(l);Vc(2);Vc(3);Vc(4);Vc(5);Vc(6) 

FOR  j%=l  TO  3 
FOR  k%=j%+l  TO  4 

IF  Vc(j%)>Vc(k%)  THEN 
x=Vc(k%) :  p=plot(k%) 

Vc(k%)=Vc(j%) :  plot(k%)=plot(j%) 

Vc(j%)=x :  plot(j%)=p 
END  IF 
NEXT  k% 

NEXT  j% 

■PRINT  "#";Vc(l);Vc(2);Vc(3);Vc(4);Vc(5);Vc(6) 
yhl=Vc(2) :  xhl=plot(2) 
yh2=Vc(3) :  xh2=plot(3) 

IF  SGN(y2-yl)oSGN(yh2-yhl)  THEN 
'  Make  sure  segment  is  going  in  same  direction  as  line  of  integration! 

SWAP  xhl,  xh2 
SWAP  yhl,  yh2 
END  IF 
END  IF 
elem%=i% 

EXIT  SUB 

nextelem:  menu0%=MENU(0) :  IF  menu0%>0  THEN  WINDOW  1  :  EXIT  SUB 
NEXT i% 
elem%=-l 
END  SUB 


m  o  a)  hb  m  m  i  cixi  9  \ 


This  program  creates  a  finite  element  model  of  two  iron  bars  (upper  and  lower)  with  an 
iron  stator  and  rotor  above  and  below.  The  dimensions  of  the  bars  and  air  spaces  and 
current  densities  are  read  in  from  a  text  file.  A  sample  text  file  is  shown  below: 

Contra  Rotating  Demo  Motor  Br=.  76  Tesla,  Ibar  =35000 


WO=. 02797 
GO=. 00127 
Hl=. 04191 
H2=. 032258 
Gl=. 00127 
G2=. 00254 


Width  of  bars 

Gap  between  bar  and  radial  symmetry 
Height  of  lower  bar 
Height  of  upper  bar 

Gap  between  Iron  core  and  lower  bar 
Gap  between  Iron  bars 
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G3=. 00127  '  Gap  between  Iron  core  and  upper  bar 

HC=. 015  ‘  Height  of  Iron  core,  top  and  bottom 

CURDl=-8764l9. 0  ’  Current  density  of  lower  bar 

CURD2=1 133616. 0  ‘  Current  density  of  upper  bar 

by=.76  '  Radial  Flux  density 

dn=699  ‘  Number  of  nodes/meter 

The  iron  bar  model  is  "hard  wired"  into  the  program.  Dimensions  and  current  densities 
can  be  changed,  but  not  the  structure  of  the  model. 

'  CREATE  INPUT  MESH  FILE  FOR  RICHARD'S  FE  PROG  7/1 8/89 
’  TWO  IRON  BARS  ON  TOP  ONE  ANOTHER  WITH  AXIAL  CURRENT 
AND  RADIAL  FLUX 
'  SIDES  OF  BARS  SEPARATED  BY  AIR 
'  WITH  SYMMETRY  +  CONSTANT  BOUNDRY  CONDITION. 

i 

'  Read  input  File,  Then  creates  or  overwrites  output  file  named  FileNameS  as  #1 

COMMON  WhoAmI$EileName$Printf%,menuO%jnenul% 

'  Who  Am  IS  is  Program  name  to  return  to  (the  shell  program) 

'  FileNameS  is  the  name  of  the  Random  acess  file  to  be  created  by 
this  program.  (The  shell  program  needs  to  know  this) 

'  Pnntf%  is  a  flag  indicating  whether  file  prints  are  to  the  screen 
or  to  a  file.  (The  shell  program  needs  to  remember  this) 

'  menuO%,  menul%  points  where  to  return  to  in  the  main  program. 

'  DO  NOT  DEFINE  OR  CHANGE  menu0%  OR  menu  1  %. 

MENU  RESET 

Begin:  CNS=nLESS(l,TEXT) 

IF  CNS=""  THEN  CHAIN  WhoAmIS 
OPEN  CN$  FOR  INPUT  AS  #3 
pathmark%=l 

WHILE  INSTR(paihmaik%,CNS,":")>0 
patlunark%=INSTR(pathmark%tCN$,":,')+ 1 
WEND 

FileName$=MID$(CN$,pathmark%.pathmark%+25)+".mesh'' 

PRINT  FileNameS 

changeCursor  4 

TWO%=0 :  TGO%=0 :  TH1%=0:  TH2%=0 :  TGl%=0  :  TG2%=0:  TG3%=0:  THC%=0 
TCURD1%=0 :  TCURD2%=0:  TBy%=0:  Tdn%=0 

WHILE  NOT(EOF(3)) 

LINE  INPUT  #3,  AS 
NC%=INSTR(A$,"='') 

IFNC%>1  THEN 
V=VAL(MIDS(AS,NC%+1)) 

AS=UCASE$(MIDS(  AS.l  J9C%- 1)) 

PRINT  AS,V 

IF  AS="WO"  THEN  WO=V  :  TW0%=1 
IF  A$=”GO"  THEN  GO=V  :  TG0%=1 
IF  A$="Hr  THEN  H1=V  :  TH1%=1 
IF  AS="H2"  THEN  H2=V  :  TH2%=1 
IF  A$="Gr  THEN  G1=V  :  TG1%=1 
IF  A$*"G2"  THEN  G2=V  :  TG2%=1 
IF  A$="G3"  THEN  G3=V  :  TG3%=1 
IF  AS="HC"  THEN  HC=V  :  THC%=1 
IF  AS*"CURDr  THEN  CURD1=V  :  TCURD1%=1 
IF  A$="CURD2"  THEN  CURD2=V  :  TCURD2%=  1 
IF  ASa'BY"  THEN  By=V  :T  By%=l 
IF  AS="DN"  THEN  dn=V  :T  dn%=l 
END  IF 
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WEND 
CLOSE  #3 

TEST9b=TWO%*TGO%*TH  1  %*TH2%*TG  1  %*TG2%*TG3%*THC%*TCURD1  %*TCURD2%*TBy% 
*Tdn% 

IF  TEST%=0  THEN  PRINT  "ERROR  IN  INPUT  FILE  Variable  Not  Defined"  :  INPUT  x 

IF  WO*GO*H1*H2*G1*G2*G3*HC=0  THEN  PRINT  "ERROR  IN  INPUT  FILE  Bad  Variable" : 

INPUT  x 

Vconst=By*(GO+WO+GO) 

PRINT  "Vconst  =";Vconst 

NGOe-CINT  (dn*GO+  .5) 

NWO=UNT(dn*WO+.5) 

NHC=CINT(dn*HC+.5) 

NGl=CINT(dn*Gl+.5) 

NHl=CINT(dn*Hl+.5) 

NG2=CINT(dn*G2+.5) 

NH2=CINT(dn*H2+.5) 

NG3=CINT  (dn*G3+.5) 

PRINT 

PRINT  NGOJ4WOJSIGO 

PRINT 

PRINT  NHC 

PRINT  NG3 

PRINT  NH2 

PRINT  NG2 

PRINT  NH1 

PRINT  NG1 

PRINT  NHC 

loop:  nl=l 
n2=nl+NHC 
^=n2+NGl 
...select  n4=n3+NHl 
n5=n4+NG2 
n6=n5+NH2 
while2:  n7--n6+NG3 
scase:  n8=n7+NHC 
n9=n8*NGO+l 
nlO=n9+NHC 
nll=nl(>fNGl 
nl2=nll+NHl 
nl3=nl2+NG2 
nl4=nl3+NH2 
n!5=nl4+NG3 
N16=nl5+NHC 
nl7=N16+n8*(NWO-l)+l 
nl8=nl7+NHC 
nl9=nl8+NGl 
n20=nl9+NHl 
n21=n20+NG2 
p22=n21+NH2 
n23=n22+NG3 
N24=n23+NHC 

n25=N24+n8*(NGO-lHl 

n26=n25+NHC 

n27=n26+NGl 

n28=n27+NH  1 

n29=n28+NG2 

n30=n29+NH2 

n31=n3CM-NG3 
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N32=n3 1+NHC 
PRINT 

PRINT  n8,Nl6.N24,N32 
Finished:  PRINT  n7  ^  1 5  ji23  ^i3 1 
PRINT  n6,nl4ji224i30 
PRINT  n5^13,n21,n29 
PRINT  n4,nl2,n20,n28 
PRINT  n3j»IU19,n27 
GoHome:  PRINT  1124)1041184126 
PRINT  nl4i94il7,n25 

n%=INSTR(CN$,":") :  NC%=0 
WHILE  n%>0 
errorl:  NC%=n% 
n%=INSTR(n%+ 1  ,CN$,": ") 

WEND 

CLOSE  #1 

OPEN  "0\#1  J?ik5Name$,8 
CLOSE  #1 

OPEN  "R",#  1  J;ileName$,8 

NAME  FileNameS  AS  FileNameS,  "MESH" 

FIELD  #1. 2  AS  aa2$,  2  AS  ab2S,  2  AS  ac2$,  2  AS  ad2$ 
FIELD  #1, 2  AS  ba2$,  2  AS  bb2$,  4  AS  bc4S 
error2:  FIELD  #1, 4  AS  ca4S,  4  AS  cb4$ 

Nodes%=N32 

PRINT  NODES  =";Nodes%;"***" 


n&=  1  '  Node  #  1  is  stored  in  Record  #2,  Node#2  in  Record#3,  and  so 

dx=GO/NGO 

x=-dx 

FOR  i%=nl  TO  n9  STEP  n8 
x=x+dx 
dy=HC/NHC 
y=-dy 

FOR  j%=nl  TO  n2 
Abort:  n&=n&+l 
y=y+dy 

LSET  ba2$=MKI$(n&-l)  ’Nu% 

LSET  bc4$=MKS$(0)  '  Vc 
PUT  #l,n& 

LSET  ca4$=MKS$(x)  ’x 
LSET cb4$=MKSS(y)  'y 
PUT  #1,  n&+Nodes% 

NEXT  j% 

dy=Gl/NGl 
FOR  j%=n2+l  TO  n3 
n&=n&+l 
y=y+dy 

LSETba2$=MKI$(n&-I)  ’Nu% 

LSET bc4$=MKS$(0)  ’Vc 
PUT  #1,  n& 

LSET  ca4$=MKS$(x)  ’x 
LSET  cb4$=MKS$(y)  ’y 
PUT  #1,  n&+Nodes% 

NEXT  j% 


dy=Hl/NHl 
FOR  j%=n3+l  TO  n4 
n&=n&+l 
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y=y+dy 

LSET  ba2S=MKIS(n&-l)  ’Nu% 
LSET  bc4$=MKSS(0)  ’  Vc 
PUT  #1,  n& 

LSET ca4$=MKS$(x)  'x 
LSET  cb4$=MKS$(y)  *y 
PUT  #1,  n&+Nodes% 

NEXTj% 

dy=G2/NG2 
F0Rj%=n4+l  T0n5 
n&=n&+l 
y=y+dy 

LSET  ba2$=MKI$(n&-l)  'Nu% 
LSET  bc4S=MKS$(0)  '  Vc 
PUT#l,n& 

LSET  ca4S=MKSS(x)  'x 
LSET  cb4$=MKSS(y)  *y 
PUT  #1,  n&+Nodes% 

NEXTj% 


dy=H2/NH2 
FOR  j%=n5+l  TO  n6 
n&=n&+l 
ysy+dy 

LSET  ba2$=MKIS(n&-l)  'Nu% 
LSET  bc4$=MKS$(0)  'Vc 
PUT  #1,  n& 

LSET  ca4$=MKS$(x)  ’x 
LSET cb4$=MKS$(y)  ’y 
PUT  #1,  n&+Nodes% 

NEXT  j% 

dy=G3/NG3 
FOR  j%=n6+l  TOn7 
n&=n&+l 
ysy+dy 

LSET  ba2$=MKI$(n&-l)  ’Nu% 
LSET  bc4$=MKSS(0)  1  Vc 
PUT  #1,  n& 

LSET ca4S=MKS$(x)  'x 
LSET  cb4$=MKSS(y)  'y 
PUT  #1 ,  n&+Nodes% 

NEXT  j% 

dy=HC/NHC 
FOR  j%=n7+l  TOn8 
n&=n&+l 
y=y+dy 

LSET  ba2$=MKIS(n&-l)  'Nu% 
LSET  bc4$=MKSS(0)  ’  Vc 
PUT  #l,n& 

LSET  ca4$=MKSS(x)  'x 
LSET  cb4S=MKS$(y)  'y 
PUT  #1,  n&+Nodes% 

NEXT  j% 

NEXT  i% 

dx=WO/NWO 

FOR  i%=n9+n8  TO  nl7  STEP  n8 
x=x+dx 
dy=HC/NHC 

y=-dy 


FOR  j%=nl  TO  n2 
n&=n&+l 
y=y+dy 

LSHT  ba2$=MKI$(n&-l)  ’Nu% 
LSET  bc4S=MKSS(0)  *Vc 
PUT  #1,  n& 

LSET  ca4$=MKS$(x)  'x 
LSET  cb4$=MKS$(y)  'y 
PUT  #1,  n&+Nodes% 

NEXTj% 

dy=Gl/NGl 
FOR  j%=n2+ 1  TO  n3 
n&=n&+l 
y=y+dy 

LSET  ba2$=MKIS{n&-l)  *Nu% 
LSET  bc4$=MKSS(0)  1  Vc 
PUT#l,n& 

LSET  ca4S=MKS$(x)  'x 
LSET  cb4$=MKSS(y)  'y 
PUT  #1,  n&+Nodes% 

NEXTj% 

dy=Hl/NHl 
FOR  j%=n3+l  TO  n4 
n&=n&+l 
yssy+dy 

LSET  ba2$=MKIS(n&-l)  'Nu% 
LSET  bc4$=MKS$(0)  'Vc 
PUT  #1,  n& 

LSET  ca4$=MKS$(x)  'x 
LSET  cb4$=MKS$(y)  'y 
PUT  #1,  n&+Nodes% 

NEXT  j% 

dy=G2/NG2 
FOR  j%=n4+l  TO  n5 
n&=n&+l 
y=y+dy 

LSET  ba2$=MKIS(n&-l)  'Nu% 
LSET  bc4S=MKSS(0)  ’Vc 
PUT#l.n& 

LSETca4S=MKSS(x)  'x 
LSET  cb4S=MKS$(y)  'y 
PUT  #1,  n&+Nodcs% 

NEXT ]% 


dy=H2/NH2 
FOR  j%=n5+l  TO  n6 
n&=n&+l 
y=y+dy 

LSET  ba2S=MKIS(n&-l)  'Nu% 
LSET  bc4S=MKSS(0)  '  Vc 
PUT  #l,n& 

LSET  ca4S=MKSS(x)  x 
LSET  cb4S=MKS$(y)  'y 
PUT  #1,  n&+Nodes% 

NEXT  j% 

dy=*G3/NG3 
FOR  j%=n6+ 1  TO  n7 
n<fe=n&+l 

L^ET  ba2S=MKIS(n&-l)  'Nu% 


LSET  bc4S=MKS$(0)  '  Vc 
PUT#l,n& 

LSET  ca4S=MKS$(x)  *x 
LSET  cb4$=MKS$(y)  'y 
PUT  #1,  n&+Nodes% 

NEXTj% 

dy=HC/NHC 
FOR  j%=n7+l  TO  n8 
n&=n&+l 
y*y+dy 

LSET  ba2$=MKI$(n&-l)  *Nu% 
LSET bc4$=MKSS(0)  ’Vc 
PUT  #1,  n& 

LSET  ca4$=MKSS(x)  x 
LSET cb4$=MKSS(y)  'y 
PUT  #1,  n&+Nodes% 

NEXTj% 

NEXT  i% 

dx=GO/NGO 

FOR  i%=n!7+n8  TO  n25-n8  STEP  n8 
x=x+dx 
dy=HC/NHC 
y*=-dy 

FOR  j%=nl  TO  n2 
n&=n&+l 

y=y+dy 

LSET  ba2$=MKIS(n&-l)  *Nu% 
LSET  bc4$=MKS$(0)  '  Vc 
PUT  #1,  n& 

LSET  ca4$=MKS$(x)  'x 
LSET cb4$=MKS$(y)  *y 
PUT  #1,  n&+Nodes% 

NEXTj% 

dy=Gl/NGl 
FOR  j%=n2+l  TO  n3 
n&=n&+ 1 
y=y+dy 

LSETba2$=MKI$(n&-l)  'Nu% 
LSET  bc4$=MKS$(0)  'Vc 
PUT  #1,  n& 

LSET  ca4S=MKS$(x)  'x 
LSET cb4S=MKSS(y)  'y 
PUT  #1,  n&+Nodes% 

NEXTj% 

dy=Hl/NHl 
FOR  j%=n3+]  TO  i>4 
n&=n&+l 
y=y+dy 

LSETba2$=MKIS(n&-l)  ’Nu% 
LSET bc4S=MKSS<0)  'Vc 
PUT  #l,n& 

LSET  ca4S=MKSS(x)  'x 
LSET  cb4S=MKSS(y)  'y 
PUT  #1,  n&+Nodes% 

NEXTj% 

dy=G2/NG2 
FOR  j%=n4+l  TOn5 
n&=n&+l 


y=y+dy 

LSET  ba2$=MK!S(n&-l)  ’Nu% 
LSET  bc4$=MKS$(0)  '  Vc 
PUT  #1,  n& 

LSET  ca4S=MKSS(x)  *x 
LSET  cb4$=MKS$(y)  'y 
PUT  #1,  n&+Nodes% 

NEXTj% 


dy=H2/NH2 
FOR  j%*n5+l  TO  n6 
n&=n&+l 
y=y+dy 

LSET  ba2S=MKIS{n&-l)  *Nu% 
LSET  bc4S=MKS$(0)  '  Vc 
PUT#l,n& 

LSET  ca4$=MKS$(x)  *x 
LSET  cb4$=MKS$(y)  'y 
PUT  #1,  n&+Nodes% 

NEXT  j% 

dy=G3/NG3 
FOR  j%=n6+l  TO  n7 
n&=n&+l 
y=y+dy 

LSETba2$=MKI$(n&-l)  ’Nu% 
LSET  bc4$=MKS$(0)  ’Vc 
PUT  #1,  n& 

LSET  ca4$=MKS$(x)  'x 
LSET  cb4$=MKS$(y)  'y 
PUT  #1,  n&+Nodes% 

NEXTj% 


dy=HC/NHC 
FOR  j%=n7+l  TOn8 
n&=n&+l 
y=y+dy 

LSETba2S=MKIS(n&-l)  'Nu% 
LSET  bc4$=MKS$(0)  Vc 
PUT  #1,  n& 

LSET  ca4$=MKS$(x)  'x 
LSET  cb4$=MKS$(y)  ’y 
PUT  #1,  n&+Nodes% 

NEXT  j% 

NEXT i% 


i%=n25 

nn&=n& 

x=x+GO/NGO 

dy=HC/NHC 

y=-dy 

FOR  j%=nl  TOn2 
n&=n&+l 
y=y+dy 

LSETba2S=MKIS(n&-nn&)  'Nu% 
LSET  bc4$=MKSS(Vconst)  'Vc 
PUT  #1,  n& 

LSET  ca4$=MKSS(x)  x 
LSET  cb4$=MKSS(y)  'y 
PUT  #1,  n&+Nodcs% 

NEXT  j% 


dy=Gl/NGl 


FOR  j%=n2-fl  TO  n3 
n&*n&+l 
y=y+dy 

LSET  ba2$=MKIS(n&-nn&)  'Nu% 
LSET  bc4$*MKSS(V  const)  'Vc 
PUT  #1,  n& 

LSET  ca4S=MKS$(x)  x 
LSET  cb4S=MKSS(y)  'y 
PUT  #1,  n&+Nodes% 

NEXTj% 

dy=Hl/NHl 
FOR  j%=n3+l  TO  n4 
n&=n&+l 
ysy-Kly 

LSET  ba2S=MKIS{n&-nn&)  ’Nuft 
LSET  bc4$=MKS$(Vconst)  *Vc 
PUT  #l,n& 

LSET  ca4S=MKS$(x)  'x 
LSET  cb4$=MKSS(y)  ’y 
PUT  #1,  n&+Nodes% 

NEXT  j% 

dy=G2/NG2 
FOR  j%=n4+l  TO  n5 
n&=n&+l 
y=y+dy 

LSET  ba2$=MKIS(n&-nn&)  'Ntrife 
LSET  bc4S=MKSS(Vconst)  ’Vc 
PUT  #1,  n& 

LSET  ca4S=MKSS(x)  'x 
LSET  cb4$=MKS$(y)  'y 
PIT  #1,  n&+Nodes% 

NEXTj% 


dy=rI2/NH2 
FOP  j%=n5+l  TO  n6 
nA=n&+l 
y=y+dy 

LSET  ba2S=MKIS(n&-nn&)  'Nu% 
LSET  bc4$=MKS$(VconsO  '  Vc 
P'.JT  #1,  n& 

LSET  ca4S=MKS$(x)  *x 
LSET  cb4$=MKSS(y)  'y 
PJT#l,n&+Nodes% 

NEXT  j% 


dp  S3/NG3 
FO*  j%=n6-t-l  TO  n7 
r&=n&+l 
y^y+dy 

LSET  ba2S*MKIS(n&-nn&)  ’Nu% 
L  SET  bc4$»MKSS{Vconst)  *Vc 
P  IT  #1,  n& 

LSET  ca4S=MKSS(x)  'x 
LSET  cb4S=MKSS(y)  ’y 
PUT  #1,  n&+Nodes% 

NEXT  j% 

dy*HC/NHC 
FOR  j%=n7+l  TO  n8 
n&=n&+l 

y=y+dy 

LSET  ba2S=MKIS(n&-nn<t)  Nu% 
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LSET  bc4S=MKS$(Vconst)  'Vc 
PUT  #1,  n& 

LSET  ca4$=MKS$(x)  'x 
LSET  cb4S=MKS$(y)  'y 
PUT  #1,  n&+Nodes% 

NEXT  j% 

Elements%=CINT((NGONWO+NGO)*(NHC+NG  1+NH 1+NG2+NH2+NG3+NHC)) 
PRINT  "•••ELEMENTS  =";Elements%;"***" 


Elements2%=2*Elements% 

Elements3%=3*Elements% 

n&=l+3*Nodes% 

NLL=nll :  NUL=nl2 :  NLR=nl9:  mate%=3  :  curd=CURDl :  u=.001  'Lower Bar 
nreg%=l 

GOSUB  writelement 

NLL=nl3  :  NUL=nl4  :  NLR=n21 :  curd=CURD2  ’Upper  Bar 
nreg%=2 

GOSUB  wntelement 

NLL=nl  :  NUL=n2  :  NLR=n25  :  curd=0  'Lower Core 
nreg%=3 

GOSUB  writelement 

NLL=n7  :  NUL=n8  :  NLR=n3 1  '  Upper  Core 
nreg%=4 

GOSUB  writelement 

NLL=n2  :  NUL=n3  :  NLR=n26:  mate%=0  :  u=l  'Lower  Air  Gap 
nreg%=5 

GOSUB  writelement 

NLL=n4  :  NUL=n5  :  NLR=n28  '  Mid  Air  Gap 
nreg%=6 

GOSUB  writelement 

NLL=n6  :  NUL=n7  :  NLR=n31  '  Upper  Air  Gap 
nreg%=7 

GOSUB  writelement 

NLL=n3  :  NUL=n4  :  NLR«nll  '  Lower  Left  Air  Gap 
nreg%=8 

GOSUB  writelement 

NLL=n5  :  NUL=n6:  NLR=nl3  '  Upper  Left  Air  Gap 
nreg%=9 

GOSUB  writelement 

NLL=nl9 :  NUL=n20 :  NLR=n27  '  Lower  Right  Air  Gap 
nreg%=10 

GOSUB  writelement 

NLL=n2l  :  NUL=n22  :  NLR=n29  '  Upper  Right  Air  Gap 
nreg%=ll 

GOSUB  writelement 

IF  n&-(l+3*Nodes%)oElements%  THEN  PRINT  "ERROR  IN  ELEMENT  GENERATION"  :  INPUT  x 

LSET  aa2S=MKIS(Nodes%) 

LSET  ab2J=MKI$(Elements%) 

LSET  ac2S=MKIS(0)  ’TestV% 

LSET  ad2S=MKIS(0)  '  TestE% 
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PUT  #1, 1 

'  Record  #1  Stores  number  of  elements  and  nodes,  and  flags  indicating 
'  that  Potential  and  Equation  data  have  not  yet  been  computed. 

CHAIN  Who AmIS 

writeiement: 

FOR  i%=NLL  TO  NLR-n8  STEP  n8 
FOR  j%*0  TO  NUL-NLL-1 
k%*i%+j% 
n&=n&+l 

LSET  aa2$=MKI$(k%)  '  Ni% 

LSET  ab2$=MKIS(k%+ 1)  '  Nj% 

LSETac2S=MKI$(k%+l+n8)  *Nk% 

Hagain:  LSET  ad2S=MKI$(kfc+n8)  *  Nl% 

PUT  #1,  n& 

'  PRINT  n&;k%;k%+l;k%+l+n8;k%+n8;mate%;curd;u 
LSET  ba2S=MKIS<mate%)  '  mate% 

LSET  bb2$*=MKIS(0)  ’  Type%  (0=split  into  best  triangles) 

LSETbc4S=MKSS(curd)  'curd 
PUT  #1,  n&+Elements% 

LSET  ca4$=MKSS(u)  '  ul=14t 

LSET  cb4$=MKSS(u)  ’  u2=  1/p 

PUT  #1,  n&+Elements2% 

LSET  aa2$=MK!S(nreg%)  ’  Region* 

LSET  ab2$=MKIS(0)  '  Blank 

LSET  ac2$=MKI$(0)  '  Blank 

LSET  ad2S=MKI$(0)  ’Blank 

PUT  #  1 ,  n&+Elemcms3% 

NEXTj% 

NEXT  i% 

RETURN 


FINP  NOPE  ANP  ELEMENT  CONNECTIONS  PROGRAM 

This  program  searches  the  geometry  part  of  the  file  to  find  out  which  elements  are 
adjacent  to  each  other  and  which  nodes  are  connected  to  any  given  node  so  that  they  will 
be  used  in  the  equation  for  that  node.  The  element  connection  data  is  needed  by  the  main 
program  to  plot  mesh  geometry,  flux  plots  and  material  outlines.  Node  connection  data  is 
needed  by  the  equation  writer  program. 

’  Find  Element  Connections  in  2D  9/15/89. 

'  Called  by  Finite  Element  Solver  Program 

COMMON  WhoAmISJFileNaine$,Printf%,menuO‘fc,menul% 

'CLEAR  JRE(-1)+FRE(0)-150(XX)&  '  Leave  150K  for  program  heap  and  stack 

WINDOW  2,FileNameS 
OPTION  BASE  0 

'  Limitation:  16000  Nodes,  12  elements  phyically  connected  to  any  given  node 
'  Conn%(x,0)  holds  number  of  elements  specified  for  x 
'  DIM  Conn%(16000,12) 

DIM  Conn%(2000,12),nelem%(4),N%<20003) 

MENU  RESET 
MENU  1,0,1, "Status" 

MENU  1 ,1 ,1, "Return  to  Main  Program  without  saving  results’ 

MENU  2,0,1," 

MENU  3,0,1,"*  Finding  Node  Connections  •" 

IF  FileName$=""  THEN  FileNameS=FILESS(  1  ."TEXTMESH") :  OPEN  "R",#lFileName$,8 
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110=000001256637# 


FIELD  #1, 2  AS  aa2$,  2  AS  ab2S,  2  AS  ac2$.  2  AS  ad2S 
FIELD  #1, 2  AS  ba2$,  2  AS  bb2$,  4  AS  bc4$ 

FIELD  #1,4  AS  ca4S,4  AS  cb4S 

GET  #1,1 

Nodcs%=CVI(aa2S) :  Elements%=CVl(ab2$) 

Test V %=C VI(ac25) :  TestE%=CVI(ad2$) 

IF  menu0%o2  AND  (TestE%=l  OR  TestE%=3)  THEN  GoHomc 

WINDOW  2,"Finding  Nodes  connected  to  Element ", (80,40)-(400, 70),  1 

'Element  Ni  Nj  Nk  Nl" 
nijkJpt&= 1  +3*Nodes% 

GET  #1,  nijklpufc 

FOR  inijkl%=l  TO  Elements**; 

PRINT  inijkl% 

GET#1 

nelem%(0)=CVI(aa2$) 
nelem%(  1  )=CVI(ab2$) 
ne!em%(2)=CVI(ac2$) 
nelem%(3)=CVI(ad2$) 

N%(inijkl%,0)=nelem%(0) 

N%(inijkl%,l)=nelem%(l) 

N%(inijkl%.2)=nelem%(2) 

N%(inijkl%3)=nelem%(3) 

FOR  j%=0  TO  3 
Ni%=nelem%(j%) 

IF  Ni%>0  THEN 
numi%=Conn%(Ni%,0) 

IF  numi%=12  THEN  crrorl 

numi%=numi%+l 

Conn%(Ni%,0)=numi% 

Conn%(Ni%,numi%)=inijkl% 

END  IF 
NEXT  j% 

IF  MENU(0)  THEN  Abort 
NEXT  inijld% 

WINDOW  2,"Finding  Elements  connected  to  Element",(80,40)-(400,70),1 
nabcdpt&=l+3*Nodes%-t4*Elements% 

GET  #l,nabcdpt& 

FOR  i%=l  TO  Elements% 

PRINT  i% 

loop:  LSET  aa2$=MKI$(0) 

LSET  ab2S=MKI$(0) 

LSET  ac2$=MKIS(0) 
mselect:  LSET  ad2S=MKIS<0) 
nelem%(0)=N%(i%,0) 
nelem%(l)=N%(i%,l) 
while2:  nelem%(2)=N%(i%3) 
scase:  nelem%(3)=N%(i%,3) 
nelem%(4)=nelem%(0) 
jj%=3 

IF  nelem%(3)=0  THEN  nelem%(3)=nelem%(0) :  ij%=2  '  Triangular  Element 
FOR  j%=0  TO  jj% 
m%=0 

Ni%=nelem%0%) 

Nj%=nelem%(j%+  1 ) 

FOR  k%=l  TO  Conn%(Ni%,0) 

Nk%=Conn%(Ni%Jc%) 

IF  Nk%oi%  THEN 
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FOR  I%=0TO3 
IF  N%(Nk%4%)=Nj%  THEN 
m%=m%+l 
IF  m%=2  THEN  aror2 
IF  j%=0  THEN  LSET  aa2$=MKI$(Nk%) 

IF  j%=l  THEN  LSET  ab2S=MKI$(Nk%) 

IFj%=2  THEN  LSET  ac2$=MKIS(Nk%) 

IF  j%=3  THEN  LSET  ad2$=MKI$(Nk%) 

END  IF 
NEXT  1% 

END  IF 
NEXT  k  % 

IF  MENU(O)  THEN  Abort 
NEXTj% 

PUT  #1 

NEXT  i% 

Finished: 

GET  #1,1 

IF TesiE%=0  THEN  LSET  ad2$=MKIS(l)  ELSE  LSET  ad2S=MKIS(3) 

PUT  #1, 1 


GoHome: 

WINDOW  CLOSE  2 
WINDOW  1 
CHAIN  Who AmlS 


enorl: 

WINDOW  CLOSE  2 
WINDOW  1 

PRINT  "Error  in  Mesh  generation!" 

PRINT  "More  that  12  elements  physically  connected  to  node" ;Ni% 
PRINT 

FOR  i%=l  TO  12 
HUNT  Conn%(Ni%  j%) 

NEXT i% 

INPUT  x 
GOTO  Abort 


enror2: 

WINDOW  CLOSE  2 
WINDOW  1 

PRINT  "Mesh  generation  error!" 

PRINT  "Segment  from";Ni%;"  to";Nj%;"  in  element";i%; 
PRINT "  is  bordered  by  two  elements:  (only  one  allowed)" 
IF  j%=0  THEN  PRINT  CVI(aa2S);"  and";Nk% 

IF  j%=l  THEN  PRINT  CVI(ab2S);"  and";Nk% 

IF  j%=2  THEN  PRINT  CVI(ac2$);"  and";Nlc% 

IFj%=3  THEN  PRINT  CVI(ad2S);"  and"^Ik% 

INPUT  x 
GOTO  Abort 


Abort 
menu0%=0 
GOTO  GoHome 


WRITE  EQUATIONS  PROGRAM 

This  program  takes  geometry,  current,  and  boundary  condition  data  from  the  file  and 
generates  the  equations  needed  for  calculate  a  solution.  These  equations  are  stored  in  the 
file  in  a  compact  form  to  be  used  by  the  solver. 
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'  Write  Equations  in  2D  9/15/89. 

'  Called  by  Finite  Element  Shell  Program 


COMMON  WhoAmI$,FileName$,Printf%,inenuO%,menul% 

‘CLEAR  ,FRE(-  1)+FRE(0)- 150000&  '  Leave  150K  for  program  heap  and  stack 

WINDOW  lfUeNameS 
OPTION  BASE  0 

DIM  v&(2000),pnt&(8000)jelem%(8000) 

DIM  Nu%(2000),Vc(2000),x(2000),y(2000) 

DIM  Ni%(20003),ul(2000) 

DIM  mate%(2000),Type%(2000),Curd(2000) 

imax%=2000  '  Maximum  number  of  nodes  SET  EQUAL  TO  DIMENSIONED  ! 

MENU  RESET 

Begin:  MENU  1,0,1, "Status" 

MENU  1,1,1,’Retum  to  Main  Program  without  saving  results" 

MENU  2,0,1," 

MENU  3,0,1,”*  Writing  Equations  *" 

IF FileNameS=""  THEN  FileName$=FILESS(l,"TEXTMESH”) :  OPEN  "R",#lrFileNameS,8 

menu0%=2 

U0=.00000l256637# 

FIELD  #  1 . 2  AS  aa2$,  2  AS  ab2$.  2  AS  ac2S ,  2  AS  ad2$ 

FIELD  #1, 2  AS  ba2$,  2  AS  bb2$,  4  AS  bc4S 
FIELD  #  1 , 4  AS  ca4$,  4  AS  cb4$ 

OPEN  "FE  .scrap"  FOR  OUTPUT  AS  #3 
CLOSE  #3 

OPEN  "FE.scrap"  AS  #3  LEN=12 

FIELD  #3, 2  AS  da2$,  2  AS  db2S,  4  AS  dc4$,  4  AS  dd4$ 

WINDOW  2,"Reading  File  Information", (80.40M400 ,70),  1 
•WINDOW  2,"Reading  File  Information",(10,40M500335),l 
GET  #1, 1 

Nodes%=CVI(aa2$) :  Elements%=CVI(ab2$) 

Test V%=C VI(ac2$) :  TeslE%=CVI(ad2$) 

IF  menu0%o2  AND  (TestE%=2  OR  TestE%=3)  THEN  GoHome 

xypt&=  1  +Nodes% 
nijklpt&=  l+3*Nodes% 
mcurdpt&=  1  +3*Nodes%+Elements  % 
typept&s  1 +3*Nodes%+Elements% 
upt&=l+3*Nodes%+2*Elements% 
equpt&=2+3*Nodes%+5*Elements% 

StartVequ&=equpt&+4  '  •**  Potential  equations  start  here 

‘  PRINT  "Reading  Nu,  Vc" 

GET  #1,1 

FOR  inuvc%=l  TO  Nodes% 

GET#1 

Nu%(inuvc%)=CVl(ba2$) 

Vc(inuvc%)=CVS(bc4S) 

IF  MENU(O)  THEN  Abort 
NEXT  inuvc% 

'  PRINT  "Reading  X,  Y" 

GET  #1,  xyptA 

FOR  ixy%=  1  TO  Nodes  % 

GET  #1 
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x(ixy%)=CVS(ca4$) 

y(ixy%)=CVS(cb4$) 

IF  MENU(O)  THEN  Abort 
NEXT  ixy% 

'  print  "Reading  Nijkl" 

GET#1,  nijldpufe 

FOR  inijkl%=  1  TO  Elements% 

GET  #1 

Ni%(inijkl%,0)=CVI(aa2$) 

Ni%(inijkl%,l)=CVl(ab2$) 

Ni%(inijkl%,2)=CVI(ac2$) 

Ni%(inijkl%3)=CVI(ad25) 

IF  MENU(O)  THEN  Abort 
NEXT  inijkl% 

loop: 

'  print  "Reading  Mate%,  Type%,  and  Curd" 

'  Assumes  xy  and  nijkl  have  been  read  in. 
mselect  GET  #1.  mcurdpufc 
FOR  imcurd%=l  TO  Elements'*? 

GET#1 

while2:  mate%(imcurd%)=CVI(ba2$) 
scase:  Type%(imcurd%)=CVI(bb2$) 

Curd(imcurd%)=C  V  S(bc4$) 

IF  Type%(itype%)=0  THEN 

nl%=Ni%(i%,0):  n2%=Ni%(i%,l) :  n3%=Ni%(i%,2) :  n4%=Ni%(i%3) 

xi=x(nl%) :  yi=y(nl%) 

xj=x(n2%) :  yj=y(n2%) 

xk=x(n3%):  yk=y(n3%) 

xl=x(n4%):  yl=y(n4%) 

d  1  =(xk-xi)A2+(yk-yi)A2 

d2=(xl-xj)A2+(yl-yj)A2 

'  Type%=l 

'  IF  dl<d2*.9999  THEN  Type%=2 
Type%=2 

IF  dl>d2*  1.0001  THEN  Type %=1 
LSET  bb2S=MKI$(Type%) 

PUT  #1,  imcurd%+mcurdpt& 

T  ype%(imcurd%)=Type% 

END  IF 

IF  MENU(O)  THEN  Abort 
NEXT  imcurd% 

'  Reading  Permeability 
GET  #1,  upt& 

FOR  i%=l  TO  Elements% 

GET  #1 

'  uO  only  used  for  linear  elements  (mate=0,l,2) 

'  ulou2  only  if  nonlinear  element  (mate>2) 
ul(i%)=CVS(ca4$) 

Finished:  u2=CVS(cb4$) 

IF  ul(i%)ou2  AND  mate%(i%)<3  THEN  PRINT  "ERROR  IN  MESH  DATA" 

IF  MENU(O)  THEN  Abort 
NEXT i% 

Go  Home: 

'  Nu%(i)=i  :  Node  is  independent 
'  Nu%{i)=j  :  VO)  =  V(i)  +  Vc(i) 

’  Nu%(i)=-n  :  NOT  YET  IMPLEMENTED  !!  V(n)  *  V(i)  +  Kn.  Kn  to  be  determined  by  solver 
(Vc(i)=Kn  ref#) 

'  Nu%0=0  :  means  node  has  a  set  potential. 

'  Note:  In  the  Solver  u{0)=l  and  v(0)=0 
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error  1: 

'  Eliminate  indirect  relationships 

'  In  the  future,  this  could  be  very  useful  for  eliminating  duplicate  nodes. 
WINDOW  2,"Investigating  Node  relationships" 

FOR  Node%=l  TO  Nodes  % 
test%=0 

j%=Nu%(Node%) :  Vc=Vc(Node%) 

WHILE  j%oNu%(j%)  AND  Nu%(j%)>0 
test%=l 
Vc=Vc+VcO%) 
j%=Nu%(j%) 

error2:  IF  MENU(O)  THEN  Abort 
WEND 

IF  test%  THEN  Nu%(Node%)=j%  :  Vc(Node%)=Vc 
'  PRINT  "combine"  ;Nu%(Node%);Vc(Node%) 

IF  MENU(O)  THEN  Abort 
NEXTNode% 

'  Link  elements  mathematically  related 

WINDOW  2,"Investigating  Element  relationships" 
vpt&=0 

FOR  i%=l  TO  Elements% 

FOR  j%=0  TO  3 

Abort  nl%=ABS(Nu%(Ni%(i%j%))) 

IF  nl%>0  THEN 
vpt&=vpt&+l 

IF  v&(nl%)=0  THEN  spt&=vpt&  ELSE  spt&=pnt&(v&(nl%)) 

pnt&(v&(nl  %))=vpt& 

v&(nl%)=vpt& 

pnt&(vptA)=spt& 

elem%(vpt&)=i% 

END  IF 
NEXT  j% 

IF  MENU(O)  THEN  Abort 
NEXT i% 


’  First  write  nodes  with  set  potential. 

nloop:  ’  These  single  line  equations  have  the  form  V(Node%)=Vc(Node%) 


WINDOW  2,"Finding  Nodes  with  set  potential" 

GET  #1,  StartVequA  '  Start  of  Node  equations 

inodes%=0  '  Number  of  equations  for  set  potential  nodes 

FOR  Node%=l  TO  Nodes%  '  First  write  nodes  with  set  potential. 

IF  Nu%(Node%)=0  THEN 
LSET  ba2$=MKIS(Node%) 

LSET  bb2$=MKIS(0)  Blank 
LSET  bc4$=MKS$(Vc(Node%)) 

PUT  #1 

inodes%=inodes%+ 1 
v&(Node%)=0 

'  PRINT  "V("+MIDS(STR$<Node%),2>+")  =’;vc(Node%) 

PRINT  Node% 

END  IF 

IF  MENU(O)  THEN  Abort 
NEXT  Node% 

StartDequ&=LOC(l>*l&  '  Start  of  Dependent  equations 
GET  #1,  Stan  VequA 
LSET  ba2$=MKIS(inodes%) 

’  LSET  bb2$=MKI$(Pnode%)  ’  Node  to  print  out  *Set  in  Main  Program* 
LSET  bc4$=MKL$(StartDequ&) 

PUT  #1,  StartVequA 
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'  Next  write  dependent  node  equations. 

’  These  single  line  equations  have  the  form  V(Node%)=V(Nu%(Node%))+Vc(Node%) 

WINDOW  2,"Finding  Dependent  Node  equations" 

GET  #1,  StartDequ&  '  Start  of  dependent  node  equations 
dnodes%=0  1  Number  of  equations  for  dependent  potential  nodes 
FOR  Node%=l  TO  Nodes% 

IF  Nu%(Node%>0  AND  Nu%(Node%)oNode%  THEN 
LSET  ba2$=MKI$(Node%) 

LSET  bb2$=MKI$(ABS(Nu%(Node%))) 

LSET  bc4$=MKSS(Vc(Node%)) 

put#i 

dnodes%=dnodes%+ 1 
v&(Node%)=0 

*  PRINT  "V("+MIDS(STR$(Node%),2)+")  =  V("+MID$(STRS(ABS(Nu%(Node%))),2>+") 
+”;vc(Node%) 

PRINT  Node% 

END  IF 

IF  MENU(O)  THEN  Abort 
NEXT  Node% 


'  Write  Independent  Node  equations 

GoHome: 

cunentnode&sl 

WINDOW  2,"Fonnulating  Equation  for  Node" 

FOR  Node%=l  TO  Nodes% 

Finished:  cl=0 :  c2=0 :  Nline%=0 
vpt&=v&(Node%) 

IF  vpt&>0  THEN 
GET  #3,  currentnodeA 
PRINT  Node% 
sptp&=0:  sptA=pntA(vpt&) 

WHILE  spt&>spq>& 

Element%=elem%(spt&) 

nl%=Ni%(Element%,0) :  n2%=Ni%(Element%,l) :  n3%=Ni%(Element%,2) : 
n4%=Ni%(Element%r3) 

Ja=Curd(Element%) 

SELECT  CASE  Type%(Element%) 

CASE  1 

nli%=nl%:  nlj%=n2%  :  nlk%=n4% 
n2i%=n2%  :  n2j%=n3%  :  n2k%=n4% 

IF  Nu%(nl%)=Node%  THEN 
testl%=l :  test2%=0  :  VcsVc^l^) 

ELSEIF  Nu%(n2%)=Node%  THEN 
testl%=2 :  test2%=l:  Vc=Vc(n2%) 

ELSEIF  Nu%(n3%)=Node%  THEN 
testl%=0:  tesl2%=2 :  Vc=Vc(n3%) 

ELSE 

testl%=3  :  test2%=3  :  Vc=Vc(n4%) 

END  IF 
CASE  2 

nli%«nl%  :  nlj%=n2%  :  nlk%=n3% 
n2i%=nl%  :  n2j%=n3%  :  n2k%=n4% 

IF  Nu%(nl  %)=Node%  THEN 
testl%=l  :  tesi2%=l :  Vc=Vc(nl%) 

ELSEIF  Nu%(n2%)=Node%  THEN 
testl%=2 :  test2%=0:  Vc=Vc(n2%) 

ELSEIF  Nu%(n3%)=Node%  THEN 
testl%=3 :  tesi2%=2:  Vc=Vc(n3%) 

ELSE 
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testl%=0:  test2%=3 :  Vc=Vc(n4%) 

END  IF 

CASEO  'Triangle-  Only  one  will  produce  nonzero  area. 
nli%=nl%  :  nlj%=n2%  :  nlk%=n3% 
n2i%=n2%  :  n2j%=n3%  :  n2k%=n4% 

IF  Nu%(nl%)=Node%  THEN 
testl%=l :  test2%=0:  Vc=Vc(nl%) 

ELSEIF  Nu%(n2%)=Node%  THEN 
testl%=2:  test2%=l :  Vc=Vc(n2%) 

ELSEIF  Nu%(n3%)=Node%  THEN 
testl%=3  :  test2%=2 :  Vc=Vc(n3%) 

ELSE 

testl%=0 :  test2%=3:  Vc=Vc(n4%) 

END  IF 
CASE  ELSE 

PRINT  "Error  in  element  types...." :  INPUT  x  :  STOP  ’  This  shouldn't  happen 
END  SELECT 
xli=x(nli%) :  yli=y(nli%) 
xlj=x(nlj%):  ylj=y(nlj%) 
xlk=x(nlk%) :  ylk=y(nlk%) 

al4=((xli-xlj)*6'l«-ylk)+{yli-ylj)*(xlk-xli))  '2* Area  positive  direction. 

'  The  direction  of  a  14  does  not  matter.  A  positive  direction  produces 
’  a  positive  Kii  if  the  element  defined  in  a  counter-clockwise  direction. 

IF  al4o0  THEN 
bli=ylj-ylk :  cli=xlk-xlj 
blj=ylk-yli :  clj=xli-xlk 
blk=yli-ylj :  clk=xlj-xli 
klii=(bli*bli+cli*cli)/al4 
Klij=(bli*blj+cli*clj)/al4 
Klik=(bli*blk+cli*clk)/al4 
k  1  jj=(b  1  j*  b  l>+c  1  j  *c  1  j)/a  1 4 
K 1  jk=(b  1  j*b  lk+c  1  j*c  1  k)/al4 
klkk=(blk*blk+clk*clk)/al4 

d=bli*bli+cli*cli+blj*blj+clj*clj+blk*blk+clk*clk 

eli=U0*Ja*al4*(b!i*bli-K:li*cli)/d 

elj=U0*Ja*al4*(blj*blj+clj*clj)/d 

elk=U0*Ja*al4*(blk*blk+clk*clk)/d 

ELSE 

klii=0:  Klij=0 :  Klik=0 
kljj=0 :  Kljk=0 
klkk=0 

eli=0:  elj=0:  elk=0 
END  IF 

IF  MENU(O)  THEN  Abort 
x2i=x(n2i%) :  y2i=y(n2i%) 

\2j=x(n2j%) :  y2j=y(n2j%) 
x2k=x(n2k%) :  y2k=y(n2k%) 

a24=((x2i-x2j)*(y2i-y2k)+(y2i-y2j)*(x2k-x2i))  ’  2*  Area  positive  direction. 

IF  a24o0  THEN 
b2i=y2j-y2k :  c2i=x2k-x2j 
b2j=y2k-y2i  r  c2j=x2i-x2k 
b2k=y2i-y2j :  c2k=x2j-x2i 
k2ii=(b2i*  b2i+c2i  *c2i)/a24 
K2ij=(b2i*b2i+c2i*c2j)/a24 
K2ik={b2i*b2k+c2i*c2k)/a24 
k2 jj=  (b2j*  b2j+c2j  *c2j)/a24 
K2jk=(b2j*b2k+c2i*c2k)/a24 
k2kk=(b2k*b2k+c2k*c2k)/a24 

d=b2i*b2i+c2i*c2i+b2j*b2j+c2j*c2j+b2k*b2k+c2k*c2k 

e2i=U0*Ja*a24*(b2i*b2i-«:2i*c2i)/d 

e2j=U0*Ja*a24*(b2j*b2j+c2j*c2j)/d 

e2k=U0*Ja*a24*(b2k*b2k+c2k*c2k)/d 

ELSE 

k2ii=0  :  K2ij=0  :  K2ik=0 
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k2jj=0 :  K2jk=0 
k2kk=0 

c2i=0 :  c2j=0 :  e2k=0 
END  IF 

kll=0:  K12=0:  k21=0:  k22=0 
IF  testl%=l  THEN 
cl=cl+eli 

kl  l=Klij :  K12=Klik 
nll%=nlj%:  nl2%=nlk% 

END  IF 

IF  tesil%=2  THEN 
cl=cl+elj 

kll=Klij :  K12=Kljk 
nll%=nli%:  nl2%=nlk% 

END  IF 

IF  testl%=3  THEN 
cl=cl+clk 

kll=Klik :  K12=Kljk 
nll%=nli%:  nl2%=nlj% 

END  IF 

IF  test2%=l  THEN 
c2=c2+e2i 

k21=K2ij :  k22=K2ik 
n21%=n2j% :  n22%=n2k% 

END  IF 

IF  iesl2%=2  THEN 
c2=c2+e2j 

k21=K2ij :  k22=K2jk 
n21%=n2i% :  n22%=n2k% 

END  IF 

IF  lesi2%=3  THEN 
c2=c2+e2k 

k21=K2ik :  k22=K2jk 
n21%=n2i% :  n22%=n2j% 

END  IF 

IF  mate%(Element%)<3  THEN 
'  PRINT  "$";Elemem%;mate%(Element%);ul(Element%) 
kl  l=kl  l*ul(Element%) 

K12=K12*ul(Elcmcnt%) 
k21=k21*ul(Element%) 
k22=k22*  u  1  (Element%) 

Element%=0 
END  IF 

IFklloO  THEN 
Nlinc%=Nlinc%+ 1 

IF  Nu%(nl  1%)=0  THEN  Vn%=nl  1  %  ELSE  Vn%=Nu%(nl  1%) 
LSET  da2$=MKIS(Vn%) 

LSET  db2$=MKI$(Element%) 

LSETdc4S=MKS$(kll) 

LSET  dd4$=MKS$(Vc(nl  1%)-Vc) 

PUT  #3 

'  PRINT  LOC(2);CVI(da2S);CVI(db2$);CVS(dc4$);CVS(dd4$) 
f-ND  IF 

IF  MENU(O)  THEN  Abort 
IF  K12o0  THEN 
Nline%=Nline%+ 1 

IF  Nu%(nl2%)=0  THEN  Vn%=nl2%  ELSE  Vn%=Nu%(nl2%) 
LSET  da2$=MKIS(Vn%) 

LSET  db2S=MKIS(Elemcnt%) 

LSET  dc4S=MKSS(K12) 

LSET  dd4$=MKSS(V  c(n  1 2%)- Vc) 

PUT  #3 
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•  PRINT  L<X(2);CVI(da2$);CVI(db2S);CVS(dc4$);CVS(dd4$) 
END  IF 

IF  MENU(O)  THEN  Abort 
IFk21oO  THEN 
Nline%=Nline%+l 

IF  Nu%(n2 1  %)=0  THEN  Vn%=n21%  ELSE  Vn%=Nu%(n21%) 
LSET  da2$=MKIS(Vn%) 

LSET  db2$=MKI$(Elemcnt%) 

LSET  dc4$=MKS$(k21) 

LSET  dd4$=MKS$(V  c(n2 1  %)- Vc) 

PUT  #3 

'  PRINT  LOC(2);CVI(da2$);CVI(db2$);CVS(dc4$);CVS((id4$) 
END  IF 

IF  MENU(O)  THEN  Abort 
IFk22oO  THEN 
Nline%=Nline%+ 1 

IF  Nu%(n22%>=0  THEN  Vn%=n22%  ELSE  Vn%=Nu%(n22%) 
LSET  da2$=MKI$(Vn%) 

LSET  db2$=MKI$(Element%) 

LSET  dc4$=MKS$(k22) 

LSET  dd4$=MKSS(Vc(n22%)-Vc) 

PUT  #3 

1  PRINT  LOC(2);CVI(da2$);CVI(db2$);CVS(dc4$);CVS(dd4$) 
END  IF 

IF  MENU(O)  THEN  Abort 


sptp&=spt& :  spt&=pnt&(spt&) 

WEND 

LSET  da2$=MKI$(Node%) 

LSET  db2$=MKI$(-l) 

LSET  dc4$=MKS$(-c  1  -c2) 

LSET  dd4$=MKSS(0) 

PUT  #3,  currentnooe& 

1  PRINT  LOC(2);CVI(da2$);CVI(db2$);CVS(dc4$);CVS(dd4$) 
currentnodeA=currentnode&+Nline%+ 1 
END  IF 
NEXT  Node% 

LSET  da2S=MKI$(-l) 

LSET  db2$=MKI$(-l) 

LSET  dc4S=MKS$(0) 

LSET  dd4$=MKS$(0) 

PUT  #3,  currentnode& 


'  Make  a  note  of  needed  tcmpory  variables 
'  Trash  info  in  v&Q.  pnt&O.  elem%0,  Nu%0,  VcQ,  CurdQ 


WINDOW  2,’Compacung  Equations" 
currentnode&=curremnode&- 1  & 
zeroS=MKS$(0) 
nsub&=0 

FOR  i%=l  TONodes% 

Vc(i%)=0 
NEXT i% 

FOR  i&=2&  TO  cunentnode& 

GET  #3,  i& 

IF  dd4$ozero$  THEN 
nsub&=nsub&+ 1  & 
pnt&(nsub&)=i& 
Vn%=CVI(da2$) 
elem%(nsub&)= V  n  % 
Vc(Vn%)=-l 
END  IF 
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IF  MENU(O)  THEN  Abort 
NEXT  i& 

'  Add  temporary  nodes  to  dependent  equations 

isub%=0 

i&=l& 

WHILE  i&<=nsub& 

GET  #3,  pnt&(i&) 
dnodes%=dnodes%+ 1 
isub%=isub%+l 
Vtest$=da2$ :  VctestS=dd4$ 

’  Write  new  dependent  equation 
LSET  ba2$=MKI$(Nodes%+isub%) 

LSET  bb2S=V  tests 
LSET  bc4S=Vctest$ 

PUT  #1,  StartDequ&+dnodes% 

'  Change  node  number 

LSET  da2$=MKI$(Nodes%+isub%) 

LSET  dd4$=zero$ 

PUT  #3,  pnt&(i&) 

'  PRINT  pnt&(i&);StartDequ&+dnodes% 
pnt£(i&)=StartDequ&+dnodes% 

'  Look  for  other  similar  nodes  to  change 
j&=i&+l 

WHILE  j&<=nsub& 

GET  #3,  pnt&(j&) 

IF  da2$=V  test$  AND  dd4$=Vctest$  THEN 
LSET  da2$=MKI$(Nodes%+isub%) 

LSETdd4$=zero$ 

PUT  #3,  pnt&(j&) 
nsub&=nsub&-l& 

FOR  k&=j&  TO  nsub& 
pn  t&(k&)=pnt&(k&+ 1  &) 
elem  %(k&)=elem%(k&+ 1  &) 

NEXT  k& 
j&=j&-l& 

END  IF 
j&=j&+l& 

WEND 

IF  MENU(O)  THEN  Abort 
i&=i&+l<& 

WEND 

StartIequ&=StartDequ&+dnodes%+l&  '  Start  of  Independent  equations 
LSET  ba2S=MKIS(dnodes%) 

LSET  bb2$=MKI$(0)  '  Not  Used  (Number  of  extra  temporary  nodes  required  in  solver) 
LSET  bc4$=MKLS(StartIequ&) 

PUT  #1,  StartDequ& 


'  Write  independent  equations  in  compact  form 
’  Trash  info  in  v&Q,  pnt&O.  elem%0,  Nu%0.  Vc 0.  CurdO 


WINDOW  2,"Wnting  Equation  for  Node" 

GET  #1,  StartlequA  1  Start  of  independent  node  equations 
inodes%s0  '  Number  of  equations  for  independent  nodes 
iUncs%=0  '  Number  of  file  lines  required  for  independent  nodes 
GET  #3, 1 

nl%=CVI(da25) :  n2%=CVI(db2S) :  K12=CVS(dc4S) 

WHILE  nl%>0 
WHILE  n2%=-l 
Node%=nl%  :  cl=K12 
GET  #3 
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nl%=CVl(da2$) :  n2%=CVl(db2$) :  K12=CVS(dc4$) 
WEND 
isub%=0 

’  PRINT  "!!!";Node%;cl 
WHILE  n2%>=0 
isub%=isub%+ 1 

Nu%(isub%)=n  1  %  :  v<fe(isub%)=n2%  :  Curd(isub%)=K12 
’  PRINT  H@";isub%jil%;n2%;K12 
GET  #3 

nl%=CVI(da2$) :  n2%=CVI(db2$) :  K12=CVS(dc4$) 
WEND 

'  Sort  out  equations  to  eliminate  duplication 
Nline%=0 

FOR  i%*l  TO  isub% 
n3%=Nu%(i%) 

IF  n3%>0  THEN 
Nline%=Nline%+l 
n4%=v&(i%) 

FOR  j%=i%+l  TO  isub% 

IF  n3%=Nu%0%)  AND  n4%=INT(v&0%))  THEN 
Nu%G%)=0 

Curd(i%)=Curd(i%)+CurdG%) 

END  IF 
NEXT  j% 

END  IF 

IF  MENU(O)  THEN  Abort 
NEXT i% 

'  Write  independent  Equation 


ilines%=ilines%+l 
inodes%=inodes%+ 1 
PRINT  inodes% 

LSET  ba2$=MKIS(Node%) 

LSET  bb2S=MKI$(Nline%) 

LSET  bc4$=MKS$(cl) 

PUT  #1,  Startlequ&+ilines% 

•  PRINT  "##";StartIequ&+ilines%,CVI(ba2$);CVI(bb2S);CVS(bc4$) 

FOR  i%=l  TO  isub% 
n3%=Nu%(i%) 

IF  n3%>0  THEN 
ilines%=ilines‘fc+ 1 
LSET  ba2S=MKIS(n3%) 

LSET  bb2$=MKIS(INT(v&(i%))) 

LSET  bc4S=MKSS(Curd(i%)) 

PUT  #1,  Startlequ&+ilines% 

1  PRINT  "##";StartIequ&+ilines%,CVI(ba2S);CVI(bb2$);CVS(bc4$) 
END  IF 

IF  MENU(O)  THEN  Abort 
NEXT i% 

'  See  if  temporary  equations  needs  to  be  writen 

IF  Vc(Node%)  THEN 
FOR  i&=l&  TO  nsubA 
IF  elem%(i&)=Node%  THEN 
ilines%=iiines%+ 1 
inodes%=inodes%+l 
GET#l,pm&(i&) 

VtestS=bb2S 
LSET  bb2S*MKI$(l) 

PUT  #1,  Startlequ&+ilines% 
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•  HUNT  •@<3)r;StartIequ&+ilines%,CVI(ba2S);CVI(bb2$);CVS(bc4S) 

ilines%=ilines%+ 1 

wri  [element:  LSET  ba2$=Vtest$ 

LSET  bb2$=MKl$(0) 

LSET  bc4$=MKS$(l !) 

PUT  #1,  Startlequ&+ilines% 

'  PRINT  "@<2>2";SlartIequ&+iIines%,CVI(ba2S);CVI(bb2$);CVS(bc4$) 
END  IF 
NEXT  i& 

IF  MENU(O)  THEN  Abort 
Hagain:  END  IF 

WEND 
CLOSE  #3 

StartBequ&=StartIequ&+ilines%+l&  '  Start  ofFlux  density  equations 
LSET  ba2$=MKI$(inodes%) 

LSET bb2S=MKIS(0)  'Blank 
LSET  bo4$=MKLS(StartBequ&) 

PUT  #1,  Startlequ&  '  inodes %,  blank%,  StartBequ& 


'  Write  Permeability  equations. 


GET#1,  StartBequA  '  Start  of  Permeability  equations. 
inodes%=0  ’  Number  of  equations  for  Permeability 
FOR  j%=l  TO  Elements% 
m%=mate%(j%) 

IF  m%>2  THEN 

nl%=Ni%(j%.0) :  n2%=Ni%(j%,l) :  n3%=Ni%(j*^) :  n4%=Ni%(j%,3) 
m  point:  SELECT  CASE  Type%(j%) 

CASE  1 

nli%=nl%  :  nlj%=n2%  :  nlk%=n4% 
n2i%=n2%  :  n2j%=n3%  :  n2k%=n4% 

CASE  2 

nli%=nl%  :  nlj%=n2%  :  nlk%=n3% 
n2i%=nl%  :  n2j%=n3%  :  n2k%=n4% 

CASE  0  'Triangle-  Only  one  will  produce  nonzero  area. 

nli%=nl%  :  nlj%=n2%  :  nlk%=n3% 
n2i%=n2%  :  n2j%=n3%  :  n2k%=n4% 

CASE  ELSE 

PRINT  "Error  in  element  types...." :  INPUT  j;  :  STOP '  This  shouldn't  happen  ! 
END  SELECT 
xli=x(nli%) :  yli=y(nli%) 
xlj=x(nlj%) :  ylj*y(nlj%) 
xlk=x(nlk%) :  ylk=y(nlk%) 

al4=((xli-xlj)*(yli-ylk>+(yli-ylj)*(xlk-xli))  '2*  Area  positive  direction. 

IF  al4oOTHEN 
bli=yij-ylk :  cli=xlk-xlj 
blj=ylk-y  li :  clj=xli-xlk 
inodes%=inodes%+ 1 

LSETaa2S=MKIS0%) 

LSETab2S=MKIS<nli%) 

LSET  ac2S=MKIS(nlj%) 

LSET  ad2S=MKIS(nlk%) 

PUT  #1 

'  PRINT  "Bn"y%;TAB(10)mli%;TAB(20);nlj%;TAB(30);nlk%;TAB(40);LOC(l) 
LSET  ca4$=MKSS<c  1  i/a  14) 

LSET  cb4S=MKSS<clj/al4) 

PUT  #1 

'  PRINT  *Bc"y%;TAB(10);CVS(ca4S);TAB(25);CVS(cb4$);TAB(40)U-OC(l) 
LSET  ca4S=MKSS(bl  i/al4) 

LSET  cb4S=MKSS(blj/al4) 

PUT  #1 
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*  PRINT  "Bb"y%;TAB(10);CVS(ca4$);TAB(25);CVS(cb4$);TAB(40);LOC(l) 

END  IF 

x2i=x(n2i%) :  y2i=y(n2i%) 
x2j=x(n2j%):  y2j=y(n2j%) 
x2k=x(n2k%) :  y2k=y(n2k%) 

a24=((x2i-x2j)*(y2i-y2k)+(y2i-y2j)*(x2k-x2i))  '  2*  Area  positive  direction. 

IF  a24o0  THEN 
b2i=y2j-y2k :  c2i=x2k-x2j 
b2j=y2k-y2i :  c2j=x2i-x2k 
inodes%=inodes%+ 1 
LSET  aa2$=MKI$(j%+Elements%) 

LSET  ab2S=MKIS(n2i%) 

LSET  ac2S=MKl$(n2j%) 

LSET  ad2S=MKIS{n2k%) 

PUT  #1 

•  PRINT  "Bn";j%+Elements%;TAB(10)^2i%;TAB(20)ui2j%;TAB(3C).n2k%;TAB(40)iOC(l) 
LSET  ca4$=MKSS(c2i/a24) 

LSET  cb4S=MKSS(c2j/a24) 

PUT  #1 

’  PRINT  "Bc"y%+Ekments%;TAB(10);CVS(ca4S);TAB(25);CVS(cb4$);TAB(40);LOC(l) 
LSET  ca4S=MKS$(b2i/a24) 

LSET  cb4$=MKSS(b2j/a24) 

PUT  #1 

’  PRINT  "Bb";j%+Eleinents%;TAB(10);CVS(ca4$);TAB(25);CVS(cb4$);TAB(40);LOC(l) 
END  IF 

IF  MENU(O)  THEN  Abort 
Idle:  END  IF 
NEXT  j% 

Endequ&=LOC(l)+l  '  End  of  Flux  density  equations 
LSET  ba2$=MKI$(inodes%) 

LSET bb2$=MKI$(0)  'Blank 
LSET  bc4$=MKL$(Endequ&) 

PUT  #1,  StartBequ& 

LSET  ca4S=MKL$(StartVequ&) 

LSET  cb4$=MKL$(StartBequ&) 

PUT  #l,equpt& 

PRINT  "Node  equations  completed",TestE% 


Finished: 

GET  #1, 1 

IF  TestE%=0  THEN  LSET  ad2S=MKIS(2)  ELSE  LSET  ad2S=MKIS(3) 
PUT  #1, 1 
CLOSE  #3 
KILL  "FE  .scrap" 

GoHome: 

WINDOW  CLOSE  2 
WINDOW  1 

readnuve:  CHAIN  Who  Am  IS 

Abort 
menuO%sO 
CLOSE  #3 
KILL  "FE  .scrap" 

GOTO  GoHome 

END 
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EQUATION  SOLVER  PROGRAM 

The  equation  solver  takes  the  equations  stored  in  the  file  and  solves  them  in  an  iterative 
fashion.  The  potential  of  each  independent  node  are  calculated  based  on  the  current 
estimated  potentials  of  the  surrounding  nodes.  This  process  is  repeated  Nodes4%  times 
and  then  the  permeability  of  each  non-linear  element  is  calculated  from  BH  lookup 
tables.  The  information  in  these  tables  is  for  1010  silicon  steel. 

There  is  no  test  for  completion.  Instead,  the  potential  of  node  Pnode%  is  printed  out  after 
every  iteration.  When  the  value  of  this  potential  stops  changing,  the  solution  is  complete 
and  the  user  manually  returns  to  the  main  program. 

'  Finite  Element  Solver  Program 
’  1  Meg  Memory  Version  9/15/89 
'  Can  compute  any  problem  with 

'  <2000  nodes,  <2000  square  elements,  and  <32500  equations 
'  Requires  only  1  Meg  computer  and  no  scratch  files 

COMMON  Who  AmI$,FileName$  ,Printf%  .men  u0%  .menu  1  % 
menu0%=0 '  Don't  return  to  Solve  command  when  returning. 

CLEAR  300000& 

’CLEAR  ,FRE(-l)+FRE(0)-50000&  '  CLEAR  Kills  COMMON  variables  ! 

PRINT  FRE(0),FRE(-  l),FRE(-2) 

WINDOW  2,"Solve  Mesh" 

DIM  v(4000),u(4000) 

DIM  vn%(4000)nl%(4000).cnst(4000) 

Begin:  DIM  v%( 18000)kij(l  8000) ,ul%(18000),ck( 18000) 

DIM  ci(4000),cj(4000).bi(4000).bj(4000) 

DIM  en%(4000),v  1  %(4000),v2%(4000),v3%(4000) 

MENU  RESET 
MENU  1.0,1, "Status" 

MENU  1,1,1, "Return  to  Main  Program" 

MENU  U,0,"-" 

MENU  13.1, "Abort  without  saving  results" 

MENU  2,0,1, ’ 

MENU  3,0,1,"*  Solving  Equations  *" 

IFFileName$=""  THEN  FileName$=FILES$(  1  ."TEXTMESH") :  OPEN  "R”,#l,FiIeName$.8 
U0=  000001256637# 

FIELD  #1. 2  AS  aa2$.  2  AS  ab2$,  2  AS  ac2S,  2  AS  ad2S 
FIELD  #1. 2  AS  ba2$,  2  AS  bb2$,  4  AS  bc4S 
FIELD  #  1 , 4  AS  ca4$,  4  AS  cb4$ 

GET  #1,  1 

Nodes%=CVI(aa2$) :  Elements%=CVI(ab2S) 

PRINT  Nodes%£lemems% 

equpt&=2+3*Nodes%+5*Elements% 

GET  #l,equpt& 

StartVequ&=CVL(ca4$) 

GET  #1,  Start  VequA 
Vnode%=CVI(ba2$) 

Pnode%=CVI(bb2S) 

StartDequ&=CVL(cb4S) 

GET  #1,  StanDcquA 
Dnode%=CVI(ba2S) 

Startiequ&=CVL(bc4S) 
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GET#l,StartIequ& 

Indnodc%=CVI(ba2$) 

StartBequ&=CVL(bc4S) 

GET  #1,  StartBequ& 

Bnodc%=CVI(ba2$) 

Endequ&=CVL(bc4$) 

B  node  1%=B  node  %-l 

PRINT  Heading  in  initial  Node  potentials  and  peimabilities" 

vpt&=l+2*Nodes% 

v(0)=l 

GET  #1,  vpt& 

FOR  i%=l  TO  Nodes% 

GET#1 

v(i%)=CVS(ca4$) 

IF  MENU(0)=1  THEN  CHAIN  WhoAmIS 
NEXT  i% 

upt&*=l+3*Nodes%+2*Elements% 

u(0>=l 

FOR  i%=l  TOElements% 

GET  #1,  i%+upt& 
u(i%)=C  V  S(ca4S) 
u(i%+Elements%)=CVS(cb4$) 
loop:  IF  MENU(0>=1  THEN  CHAIN  WhoAmIS 
NEXT  i% 

mselecL  PRINT  "Reading  Set  Node  Potentials" 

GET  #1,  Start  VequA 
PRINT  StartVequ&,Vnode% 
while2:  FOR  i%=l  TO  Vnode% 
scase:  GET  #1 
v(CVI(ba2S))=CVS(bc4S) 

IF  MENU(0)=1  THEN  CHAIN  WhoAmIS 
'  PRINT  "V"+MID$(STR$(CVI(ba2$)),2}+"=";CVS(bc4$) 
NEXT i% 

'  Reading  dependent  nodes 

This  is  necessary  because  if  potentials  are  already  defined,  then 
'  then  some  other  temporary  potentials  must  be  defined. 

GET  #1 ,  StartDequ& 

FOR  i%=l  TO  Dnode% 

GET  #1 

v(CVI(ba2S))=v(CVI(bb2S)>+CVS(bc4S) 

IF  MENU(1)=3  THEN  CHAIN  WhoAmIS 
NEXT  i% 

PRINT  "Number  Crunching" 

Nodes4%=SQR(Nodes%y3+.5 


npos%=l 

ipos&=-l 

GET  #1,  StartlequA 
FOR  i%«l  TO  Indnode% 
GET  #1 

npos%=npos%+l 

Vnum%=€VI(ba2$) 

nline%*CVI(bb2$) 

Finished:  Num=CVS(bc4S) 

vn%(npos%)=Vnum% 

nI%(npos%)*nline% 

cnst(npos%)=Num 
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Den=0 

GoHome:  FOR  j%=l  TO  niine% 
GET#1 

ipos&=ipos&+l 
v%=CVI(ba2$) 
u%s=CVl(bb2$) 
error  1:  ck=CVS(bc4$) 
kij=u(u%)*ck 
v%(ipos&>v% 
ul%(ipos&)=u% 
ck(ipos&)=dc 
kij(ipos&)=kij 
Den=Den+kij 
Num=Num+v(v%)*kij 
NEXT  j% 

IF  MENU(0)=  1  THEN  GoHome 
v(Vnum%)=Num/Den 
NEXT  i% 

error2:  iequ&=ipos& 


GET#1,  StartDequA  '  dependent  potentials. 
FOR  i%=l  TO  Dnode% 

GET#1 

npos%=npos%+l 

vn%(npos%)=CVI(ba2$) 

nl%(npos%>=CVI(bb2$) 

cnst(npos%)=CVS(bc4$) 

IF  MENU(0)=  1  THEN  GoHome 
NEXT i% 

GET  #1,  StartBequ& 

Abort;  FOR  bpos%=0  TO  Bnodel% 

GET  #1 

en%(bpos%)=CVI(aa2$) 
v  1  %(bpos%)=CVI(ab2S) 
v2%(bpos%)=CVI(ac2$) 
v3%(bpos%)=CVI(ad2$) 

GET  #1 

ci(bpos%)=CVS(ca4$) 

cj(bpos9fc)=CVS(cb4$) 

GET  #1 

bi(bpos%)=CVS(ca4$) 

bj(bpos%)=CVS(cb4$) 

IF  MENU(0)=1  THEN  GoHome 
NEXT  bpos% 


nloop: 

FOR  Eioop%-l  TO  Nodes4%  '  Loop  throught  independent  equations  this  many  times  before  updating  u 
npos%=- 1 
ipos<t=-l 

FOR  i%=l  TO  Indnode% 
npos%=npos%+l 
V  num%= vn%(npos%) 
nline%=nl%(npos%) 

Num=cnst(npos%) 

Dcn*0 

FOR  j%=1  TOnUne% 
iposA=ipos&+l 
kij=kij(ipos&) 

Den=Don+kij 

Num=Num+v(v%(tpos&))*kij 
NEXT  j% 

IF  MENU(0)*1  THEN  GoHome 
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v(Vnum%)=Num/Den 
NEXT i% 

PRINT  v(Pnode%) 

NEXT  Eloop% 

'  PRINT  v(Pnode%);"  Time  =”;TIMER-xun& 

'  B  equations  may  require  these  dependent  potentials. 

FOR  i%=l  TO  Dnode% 
npos%=npos%+l 

v(vn%(npos%))=v(nl%(npos%)}+cnst(npos%) 

IF  MENU(0)=1  THEN  GoHome 

NEXT  i% 

'  Given  BA2  find  u=l/p. 

FOR  bpos%=0  TO  Bnodel% 
e%=en%(bpos%) 
u=u(e%) 

v3=v(v3%(bpos%)) 
v  13=v(v  1  %(bpos%))-v3 
v23=v(v2%(bpos%))-v3 
bx=vl3*ci(bpos%)+v23*cj(bpos%) 
by=v  13*bi(bpos%)+v23*bj(bpos%) 

FINDUB  (bx*bx+by*by)*u,u(e%) 

IF  MENU(0)=1  THEN  GoHome 

NEXT  bpos% 

FOR  i&=0  TO  iequ& 
kij(i&)=u(ul  %(i&))*ck(i&) 

NEXT  i& 

GOTO  nloop 


GoHome: '  STOP 

PRINT  "Exiting  early,  saving  results  so  far..." 

IF  MENU(  1  )=3  THEN  CHAIN  WhoAmIS  ’  This  person  really  wants  to  leave  ! 
Finished: 

'  Writing  dependent  nodes 
GET  #1,  StartDequ& 

FOR  i%=l  TO  Dnode% 

GET  #1 

v(CVI(ba2$)>=v(CVI{bb2S)>+CVS(bc4S) 

IF  MENU(1>=3  THEN  CHAIN  WhoAmIS  ’  This  person  really  wants  to  leave ! 
NEXT  i% 

vpt&= 1  +2*Nodes% 
upt4=l+3*Nodes%+2*Elements% 

FOR  i%=l  TO  Nodes% 

LSET  ca4$=MKSS(v(i%)) 

PUT  #1,  i%+vpt& 

IF  MENU(1)=3  THEN  CHAIN  WhoAmIS  ’  This  person  really  wants  to  leave  1 
NEXT i% 

FOR  i%=l  TO  Elements% 

LSET  ca4$=MKSS(u(i%)) 

LSET  cb4$aMKSS(u(i%-«-Elements%)) 

PUT  #1,  i%+upt& 

IF  MENU(1)=3  THEN  CHAIN  WhoAmIS  '  This  person  really  wants  to  leave  ! 
NEXT i% 

MENU  RESET 
CHAIN  WhoAmIS 
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•  Given  H*^iO*B=B*B/nr  Find  H*nO/B=14ir 

SUB  FINDUB(x.Y)  STATIC 
IF  x<  .143337936992808#  THEN 
IF  x<  .0132912500793557#  THEN 
IF  x<  1 .50075074867 1 93D-03  THEN 
IF  x<  4.435 1949508249 ID-04  THEN 
IF  x<  1.98247562498695D-04  THEN 
M#=0 :  c#=.00025 
ELSE 

M#=  2364707969374813# :  c#=  2.03120240903375D-04 
END  IF 

FT  SF 

IF  x<  8.8 1999875548623  ID-04  THEN 
M#=  .3238463873585482# :  c#=  1.64367463 157289D-04 
ELSE 

M#=  350707216206927#  :  c#=  1. 406762 154558863D-04 
END  IF 
END  IF 
ELSE 

IF  x<  .0036096004864159#  THEN 
IF  x<  224626054 9%972D-03  THEN 
M#=  .359442139649298#  :  c#=  1 27567272560 156D-04 
ELSE 

M#=  .348432702993648#  :  c#=  1 32297335797 132D-04 
END  IF 
ELSE 

IF  x<  5.93505072258168D-03  THEN 
M#=  .331118707005512#  :  c#=  2.147939441377122D-04 
ELSE 

IF  x<  9.363599716406721D-03  THEN 
M#=  .309168583668735#  :  c#=  3.45069039508408D-04 
ELSE 

M#=  2800656938351602# :  c#=  6.175768505006809D-04 
END  IF 
END  IF 
END  IF 
END  IF 
ELSE 

IF  x<  .0531999989162207#  THEN 
IF  x<  .0246077741219678#  THEN 
IF  x<  .018370800574581#  THEN 
M#=  261834221882022#  :  c#=  8398959035446059D-04 
ELSE 

M#=  .2437078744305423# :  c#=  1.1 9289 141 772 1294D-03 
END  IF 
ELSE 

IF  x<  .0324900008730129#  THEN 
M#=  229630605100122# :  c#=  1.539301681658383D-03 
ELSE 

IF  x<  .0422077506956306#  THEN 
M#=  216099403610992#  :  c#=  1.9789304 29853 15D-03 
ELSE 

M#=  200140999544927# :  c#=  2.652498770 173752D-03 
END  IF 
END  IF 
END  IF 
FLSE 

IF  x<  .0890819988658082#  THEN 
IF  x<  6.892 10001 70673 12D-02  THEN 
M#=  .197188478098947#  :  c#=  2.80957290789999D-03 
ELSE 
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M#=  .18848271 1579081#  :  c#=  3.409583043701 52D-03 
END  IF 
ELSE 

IF  x<  .102052912146259#  THEN 
M#=  .185025762499145# :  c#=  3.71 753497771 9552D-03 
ELSE 

IF  x<  .1183359990843981#  THEN 
M#=  .184243132182679#  :  c#=  3.79740468064885D-03 
ELSE 

M#=  .1879854316924574#  :  c#=  3.35455592928614D-03 
END  IF 
END  IF 
END  IF 
END  IF 
END  IF 
ELSE 

IF  x<  2.03033999894946#  THEN 
IF  x<  .4995462478981502#  THEN 
IF  x<  281247075211 1964#  THEN 
IF  x<  .1824680002758821#  THEN 
M#=  .18911291072243# :  c#=  3.192945411 1272 14D-03 
ELSE 

M#=  .180756681252457#  :  c#=  4.71768989235%34D-03 
END  IF 
ELSE 

IF  x<  .346509840842061 1#  THEN 
M#=  .1701129239472231#  :  c#=  .0077112155037144# 

ELSE 

IF  x<  .409378303077775#  THEN 
M#=  .1632181542895672#  :  c#=  1.010032104043142D-02 
ELSE 

M#=  .155287222941354#  :  c#=  1.334707225758943D-02 
END  IF 
END  IF 
END  IF 
ELSE 

IF  x<  .886261429992511#  THEN 
IF  x<  .639882240674094#  THEN 
M#=  .143729663076772#  :  c#=  .0191206079227995# 

ELSE 

M#=  .128816378966854#:  c#=  2.86633535748632 1 D-02 
END  IF 
ELSE 

IF  x<  1.09471035568026#  THEN 
M#=  .1 14274469452066#  :  c#=  4.1551287096261 12D-02 
ELSE 

IF  x<  1.42663184817263#  THEN 
M#=  .100388021630704#  :  c#=  5.675292532991951D-02 
ELSE 

M#=  .0828135878631228#  :  c#=  .081825172256351# 

END  IF 
END  IF 
END  IF 
END  IF 
ELSE 

IF  x<  20.34594676792601#  THEN 
IF  x<  5.083272842696882#  THEN 
IF  x<  3.4282830007005#  THEN 
M#=  .0596364952974854#  :  c#=  .128882550351719# 

ELSE 

M#=  4.039802642433373D-02  :  c#=  .1948374661490501# 
END  IF 
ELSE 

IF  x<  9. 14977047%96492#  THEN 
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M#=  2.45323448863838 ID-02  :  c#=  .2754870542417882# 

FI  £F 

IF  x<  12.8705572079986#  THEN 
M#=  .0149570400771303#  :  c#=  .363098895519592# 
ELSE 

M#=. 00929988839034 12#:  c#=  .435909589938737# 
END  IF 
END  IF 
END  IF 
ELSE 

IF  x<  40.0624352353194#  THEN 
IF  x<  27.465932977501 1#  THEN 
M#=  5.8239747297403 ID-03  :  c#=  ^066303442472294# 
ELSE 

M#=  3.78864013587588D-03  :  c#=  .5625327077890993# 
END  IF 
ELSE 

IF  x<  121.4038573017713#  THEN 
M#=  1.34537288523743D-03  :  c#=  .660415943780379# 
ELSE 

IF  x<  9999.991009271052#  THEN 
M#=  1 .784 1 5897595796D-05  :  c#=  .821583363740207# 
ELSE 

M#=0:  c#=l 
END  IF 
END  IF 
END  IF 
END  IF 
END  IF 
END  IF 
Y=x*M#+c# 

END  SUB 
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